UtahLake_TributaryData

Details

Load data

Data were loaded from spreadsheets obtained by Utah DWQ:

# tributary data
trib <- read.csv("./Data/Raw/ul_Tribdata_wqp_processed_2021-03-02.csv")
# facility data
fac <- read.csv("./Data/Raw/ul_Facilitydata_wqp_processed.csv")
# watershed spatial data
wshds <- read_sf("./Data/GIS/UL_subwshedsWGS84_v2.shp")
# sampling activity logs
activitylog <- read.csv("./Data/Raw/UL_Trib_SampleActivities_exported 12-17-20.csv")
# groundwater data
gw <- read.csv("./Data/Raw/ul_wqp_gwdata_processed_2021-01-07.csv")
# precipitation and evaporation data
p_et <- read.xlsx("./Data/Raw/ULEFDCFlowBoundaries_WY2006-2018.xlsx", 
                  sheet = "Precip-ET", startRow = 2)
# DMR data
dmr <- read.xlsx("./Data/Raw/DMR_Analysis_v3_updated 12-7-20.xlsx", 
                 sheet = "DMR.Formatted")
dmr_LUT <- read.xlsx("./Data/Raw/DMR_Analysis_v3_updated 12-7-20.xlsx", 
                     sheet = "LUT_SiteList", startRow = 2, cols = 1:5)
dmr_Payson <- read.xlsx("./Data/Raw/dmr_payson.xlsx",
                        sheet = "DMR - 30 day avg nutrients", startRow = 3)

dmr_updated <- read.xlsx("./Data/Raw/DMR_Analysis_v4_updated 12-9-2021.xlsx",
                         sheet = "DMR.Formatted")

# Jordan River outflow hydrology
jordan_outflow <- read.xlsx("./Data/Raw/Jordan_River_Outflow.xlsx", 
                  sheet = "DailyValues20219216", startRow = 5)

# Utah lake storage capacity curve
ul_level <- read.xlsx("./Data/Raw/Copy of ULDB_Stage_Vol_SA_V3(sd-JLM).xlsx", 
                  startRow = 11, cols = c(7:10))

# 2019 and 2020 WFWQC flow data
WFWQC_flow <- read.csv("./Data/Raw/WFWQC_Flow.csv")

Data processing

Watersheds

Watersheds were divided into monitored and unmonitored and arranged in clockwise order.

wshds_monitored <- wshds %>%
  filter(Monitored_ == "YES") %>%
  # Dry Creek - Saratoga appears twice. Remove the instance with very very small area. 
  filter(OBJECTID_1 != 74)

wshds_monitored$WS_Name <- 
  factor(wshds_monitored$WS_Name, 
         levels = c("Tickville Wash", "Dry Creek - Saratoga", "Lehi Spring Creek", 
                    "American Fork River", "Timp SSD", "Lindon Drain", 
                    "Powell Slough Major", "Provo River", "Mill Race", 
                    "Spring Creek - Springville", "Hobble Creek", 
                    "Dry Creek - Spanish Fork", "Spanish Fork River", 
                    "4000 South Drain Spanish Fork", "Benjamin Slough", "Currant Creek"))
levels(wshds_monitored$WS_Name)

##  [1] "Tickville Wash"                "Dry Creek - Saratoga"         
##  [3] "Lehi Spring Creek"             "American Fork River"          
##  [5] "Timp SSD"                      "Lindon Drain"                 
##  [7] "Powell Slough Major"           "Provo River"                  
##  [9] "Mill Race"                     "Spring Creek - Springville"   
## [11] "Hobble Creek"                  "Dry Creek - Spanish Fork"     
## [13] "Spanish Fork River"            "4000 South Drain Spanish Fork"
## [15] "Benjamin Slough"               "Currant Creek"

Tributary data

Tributary dataset was filtered for:

Variables pertaining to C, N, P, or flow (discharge)
Years 2010-2020

Flow data were converted to cubic feet/sec (cfs) when units originally appeared as gallon/min, cubic meters/sec, and million gallons/day.

trib$ActivityStartDate <- as.Date(trib$ActivityStartDate)
# unique(trib$CharacteristicName)
trib$ResultMeasureValue <- as.numeric(as.character(trib$ResultMeasureValue))

## Warning: NAs introduced by coercion

trib_wq <- trib %>%
  # Set nondetects (qualifier code U) as 1/2 MDL. 
  # Qualifier code J is above detection but below reporting. Values retained as-is.
  # REVISIT: What does qualifier code Q mean? One sample for ammonium at 0.062
  mutate(ResultMeasureValue = case_when(MeasureQualifierCode == "U" ~ MethodDetectionLevel/2, 
                                        TRUE ~ ResultMeasureValue)) %>%
  filter(CharacteristicName %in% 
           c("Ammonia-nitrogen as N", "Carbon", "Nitrate as N", "Organic carbon",
             "Ammonium", "Carbonate", "Inorganic carbon", "Nitrate", "Nitrogen", 
             "Phosphorus", "Ammonia and ammonium", 
             "Inorganic nitrogen (nitrate and nitrite)", "Nitrite", "Orthophosphate", 
             "Kjeldahl nitrogen", "Nitrogen, mixed forms (NH3), (NH4), organic, (NO2) and (NO3)",
             "Organic Nitrogen", "Ammonia-nitrogen", "Phosphate-phosphorus",
             "Total Kjeldahl nitrogen")) %>%
  filter(ResultMeasure.MeasureUnitCode %in% c("ueq/L", "mg/l", "mg/l as N", "<NA>", "mg/l as P")) %>%
  select(ActivityStartDate:ResultSampleFractionText,
         OrganizationIdentifier:MonitoringLocationTypeName, LatitudeMeasure, 
         LongitudeMeasure, ResultMeasureValue, ResultMeasure.MeasureUnitCode, 
         ResultStatusIdentifier) %>%
  mutate(Month = month(ActivityStartDate), 
         Year = year(ActivityStartDate)) %>%
  filter(Year >= 2010) %>%
  droplevels() %>%
  filter(!is.na(ResultMeasureValue)) %>%
  select(-SampleDepthValue, -RelativeDepth, -ResultMeasure.MeasureUnitCode) %>%
  unite("Variable", CharacteristicName:ResultSampleFractionText) 

trib_flow <- trib %>%
  filter(CharacteristicName %in% c("Flow", "Stream flow, instantaneous", "Stream flow, mean. daily")) %>%
  droplevels() %>%
  #filter(ResultStatusIdentifier != "Provisional") %>% # can filter out provisional data
  select(ActivityStartDate:ResultSampleFractionText,
         OrganizationIdentifier:MonitoringLocationTypeName, LatitudeMeasure,
         LongitudeMeasure, ResultMeasureValue, ResultMeasure.MeasureUnitCode,
         ResultStatusIdentifier) %>%
  select(-SampleDepthValue, -RelativeDepth, -ResultSampleFractionText,
         -MonitoringLocationTypeName, -CharacteristicName) %>%
  group_by(ActivityStartDate, MonitoringLocationIdentifier, OrganizationIdentifier,
           OrganizationFormalName, MonitoringLocationName, LatitudeMeasure,
           LongitudeMeasure, ResultMeasure.MeasureUnitCode, ResultStatusIdentifier) %>%
  summarise(ResultMeasureValue = mean(ResultMeasureValue)) %>%
  mutate(Month = month(ActivityStartDate),
         Year = year(ActivityStartDate)) %>%
  filter(Year >= 2010) %>%
  droplevels() %>%
  pivot_wider(names_from = ResultMeasure.MeasureUnitCode, values_from = ResultMeasureValue) %>%
  # USGS data has both cfs and cms, always duplicated. remove cfs. 
  select(-'ft3/s')

## `summarise()` has grouped output by 'ActivityStartDate', 'MonitoringLocationIdentifier', 'OrganizationIdentifier', 'OrganizationFormalName', 'MonitoringLocationName', 'LatitudeMeasure', 'LongitudeMeasure', 'ResultMeasure.MeasureUnitCode'. You can override using the `.groups` argument.

colnames(trib_flow)[11:14] <- c("m3.sec", "cfs", "mgd", "g.min")

# convert everything to cfs
trib_flow$mgd <- trib_flow$mgd/0.646317
trib_flow$g.min <- trib_flow$g.min*0.002228017
trib_flow$m3.sec <- trib_flow$m3.sec*35.3147

trib_flow_long <- trib_flow %>%
  pivot_longer(cols = c("cfs", "mgd", "g.min", "m3.sec"), 
               names_to = "orig.units", values_to = "Flow.cfs") %>%
  drop_na(Flow.cfs)

# conversions seem to match up
ggplot(trib_flow_long, aes(x = ActivityStartDate, y = Flow.cfs, color = orig.units)) +
  geom_point(alpha = 0.5) +
  scale_y_log10() +
  scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.8)

## Warning: Transformation introduced infinite values in continuous y-axis

WFWQC_flow$ActivityStartDate <- as.Date(WFWQC_flow$ActivityStartDate)

WFWQC_flow$ResultMeasureValue[WFWQC_flow$ResultMeasureValue == "-0.153"] <- NA
WFWQC_flow$ResultMeasureValue[WFWQC_flow$ResultMeasureValue == "too deep"] <- NA
WFWQC_flow$ResultMeasureValue[WFWQC_flow$ResultMeasureValue == "no flow"] <- 0
WFWQC_flow$ResultMeasureValue <- as.numeric(WFWQC_flow$ResultMeasureValue)

flow_sites_names <- trib_wq %>%
  ungroup() %>%
  select(MonitoringLocationIdentifier, MonitoringLocationName, 
         LatitudeMeasure, LongitudeMeasure,)

WFWQC_flow_processed <- WFWQC_flow %>%
  mutate(Month = month(ActivityStartDate),
         Year = year(ActivityStartDate), 
         orig.units = "m3.sec",
         Flow.cfs = ResultMeasureValue*35.3147) %>%
  left_join(., flow_sites_names) %>%
  select(ActivityStartDate, MonitoringLocationIdentifier, OrganizationIdentifier, 
         OrganizationFormalName, MonitoringLocationName, LatitudeMeasure, 
         LongitudeMeasure, ResultStatusIdentifier, Month, Year,
         orig.units, Flow.cfs) %>%
  distinct() %>%
  drop_na(LatitudeMeasure, LongitudeMeasure)

## Joining, by = "MonitoringLocationIdentifier"

trib_flow_long <- rbind(trib_flow_long, WFWQC_flow_processed)

Find non-sampled dates

A sample type of “11” indicates a site was not sampled. A sample type of “10” indicates there was no flow. The latter should appear in the flow data as a zero; this was verified by inspecting the two datasets.

Locations and dates when sites were not sampled:

##  [1] "LOWER S FK PROVO R AT GAGING STATIION"           
##  [2] "N FK PROVO R AB CNFL / PROVO R AT WILDWOOD"      
##  [3] "CURRANT CK AB MONA RES"                          
##  [4] "Timpanogos Effluent below constructed duck ponds"
##  [5] "Dry Ck Near Utah Lake-WLA"                       
##  [6] "DRY CK EAST TRIBUTARY"                           
##  [7] "AMERICAN FK CK 2.5MI S OF AM FK CITY"            
##  [8] "Powell Slough WMA North Outfall to Utah Lake"    
##  [9] "Powell Slough WMA South Outfall to Utah Lake"    
## [10] "Saratoga Springs at Cedar Valley"

##  [1] "2011-01-26" "2014-02-19" "2012-01-19" "2012-02-15" "2012-03-21"
##  [6] "2012-04-10" "2012-06-20" "2012-07-25" "2011-10-27" "2011-11-16"
## [11] "2011-12-06" "2011-12-07" "2017-05-17" "2017-07-10" "2017-12-20"
## [16] "2017-12-22" "2018-01-10" "2018-02-15" "2018-03-13" "2018-04-19"
## [21] "2018-05-17" "2018-11-20" "2019-05-14" "2019-06-13" "2019-06-17"
## [26] "2019-07-08" "2019-08-12" "2019-09-24" "2019-09-23" "2019-10-08"
## [31] "2019-11-04" "2019-12-10"

Facility data

Tributary dataset was filtered for:

Variables pertaining to C, N, P, or flow (discharge)
Years 2010-2020

Flow data were converted to cubic feet/sec (cfs) when units originally appeared as gallon/min, cubic meters/sec, and million gallons/day.

fac$datetime <- as.POSIXct(fac$datetime)
#unique(fac$CharacteristicName)

fac_wq <- fac %>%
  # Set nondetects (qualifier code U) as 1/2 MDL. 
  # Qualifier code J is above detection but below reporting. Values retained as-is.
  # REVISIT: What does qualifier code Q mean? One sample for ammonium at 0.062
  mutate(ResultMeasureValue = case_when(MeasureQualifierCode == "U" ~ MethodDetectionLevel/2, 
                                        TRUE ~ ResultMeasureValue)) %>%
  filter(CharacteristicName %in% 
           c("Ammonia-nitrogen as N", "Carbon", "Nitrate as N", "Organic carbon",
             "Ammonium", "Carbonate", "Inorganic carbon", "Nitrate", "Nitrogen", 
             "Phosphorus", "Ammonia and ammonium", 
             "Inorganic nitrogen (nitrate and nitrite)", "Nitrite", "Orthophosphate", 
             "Kjeldahl nitrogen", "Nitrogen, mixed forms (NH3), (NH4), organic, (NO2) and (NO3)",
             "Organic Nitrogen", "Ammonia-nitrogen", "Phosphate-phosphorus",
             "Total Kjeldahl nitrogen")) %>%
  filter(ResultMeasure.MeasureUnitCode %in% c("ueq/L", "mg/l", "mg/l as N", "<NA>", "mg/l as P")) %>%
  select(MonitoringLocationIdentifier:LongitudeMeasure, ResultMeasureValue, 
         ResultMeasure.MeasureUnitCode, ResultStatusIdentifier) %>%
  mutate(Month = month(datetime), 
         Year = year(datetime)) %>%
  filter(Year >= 2010) %>%
  droplevels() %>%
  filter(!is.na(ResultMeasureValue)) %>%
  select(-SampleDepthValue, -RelativeDepth, -ResultMeasure.MeasureUnitCode) %>%
  unite("Variable", CharacteristicName:ResultSampleFractionText) 

# match date field with trib
colnames(fac_wq)[2] <- "ActivityStartDate"


fac_flow <- fac %>%
  filter(CharacteristicName == "Flow") %>%
           #c("Height, gage", "Stream flow, instantaneous", "Stream flow, mean. daily")) %>%
  droplevels() %>%
  #filter(ResultStatusIdentifier != "Provisional") %>% # can filter out provisional data
  select(MonitoringLocationIdentifier:LongitudeMeasure, ResultMeasureValue, 
         ResultMeasure.MeasureUnitCode, ResultStatusIdentifier) %>%
  select(-SampleDepthValue, -RelativeDepth, -ResultSampleFractionText, 
         -MonitoringLocationTypeName, -CharacteristicName) %>%
  group_by(MonitoringLocationIdentifier, datetime, OrganizationIdentifier, 
           OrganizationFormalName, MonitoringLocationName, LatitudeMeasure, 
           LongitudeMeasure, ResultMeasure.MeasureUnitCode, ResultStatusIdentifier) %>%
  summarise(ResultMeasureValue = mean(ResultMeasureValue)) %>%
  mutate(Month = month(datetime), 
         Year = year(datetime)) %>%
  filter(Year >= 2010) %>%
  droplevels() %>%
  pivot_wider(names_from = ResultMeasure.MeasureUnitCode, values_from = ResultMeasureValue)

## `summarise()` has grouped output by 'MonitoringLocationIdentifier', 'datetime', 'OrganizationIdentifier', 'OrganizationFormalName', 'MonitoringLocationName', 'LatitudeMeasure', 'LongitudeMeasure', 'ResultMeasure.MeasureUnitCode'. You can override using the `.groups` argument.

fac_flow$datetime <- as.Date(fac_flow$datetime)
colnames(fac_flow)[2] <- "ActivityStartDate"

colnames(fac_flow)[11:14] <- c("mgd", "cfs", "g.min", "m3.sec")

# convert everything to cfs
fac_flow$mgd <- fac_flow$mgd/0.646317
fac_flow$g.min <- fac_flow$g.min*0.002228017
fac_flow$m3.sec <- fac_flow$m3.sec*35.3147

fac_flow_long <- fac_flow %>%
  pivot_longer(cols = c("cfs", "mgd", "g.min", "m3.sec"), 
               names_to = "orig.units", values_to = "Flow.cfs") %>%
  drop_na(Flow.cfs)

Combine trib and fac data

Tributary and facility water quality and flow datasets were combined.

USGS gage data for missing discharge

USGS_flow <- readNWISdv(siteNumbers = c("10153100", "10163000"),
                         parameterCd = "00060", 
                         startDate = "2010-01-01")

USGS_flow_processed <- USGS_flow %>%
  mutate(WS_Name = case_when(site_no == "10153100" ~ "Hobble Creek",
                             site_no == "10163000" ~ "Provo River")) %>%
  select(WS_Name, Date, X_00060_00003)

colnames(USGS_flow_processed) <- c("WS_Name", "ActivityStartDate", "Flow.USGS.cfs")

Create spatial datasets for tributary and facility datasets

Intersect stations with watersheds
Join water quality data with intersected dataset
If a station is not associated with a watershed, assign it
Determine downstream stations in watersheds

Sites located downstream of diversion

UTAHDWQ_WQX-4995465
WFWQC_UT-4995467
UTAHDWQ_WQX-4995312
UTAHDWQ_WQX-4995578
WFWQC_UT-4995575
UTAHDWQ_WQX-4995578
WFWQC_UT-4995575

Sites that may be backwatered at times

WFWQC_UT-4917712
WFWQC_UT-4994958
UTAHDWQ_WQX-4995465
WFWQC_UT-4995467
WFWQC_UT-4996536
UTAHDWQ_WQX-4995210
UTAHDWQ_WQX-4995230
WFWQC_UT-4995210
WFWQC_UT-4995575

Intermittent sites

UTAHDWQ_WQX-4994804
UTAHDWQ_WQX-4994792

Effluent-dominated sites

UTAHDWQ_WQX-4996040
WFWQC_UT-4996040

Sites included in ULWQS SAP

UTAHDWQ_WQX-5919910
UTAHDWQ_WQX-4994960
UTAHDWQ_WQX-4995465
UTAHDWQ_WQX-4994804
UTAHDWQ_WQX-4996040
UTAHDWQ_WQX-4996100
UTAHDWQ_WQX-4994792
UTAHDWQ_WQX-4994950
UTAHDWQ_WQX-4995120
UTAHDWQ_WQX-4996540
UTAHDWQ_WQX-4995210
UTAHDWQ_WQX-4995230
UTAHDWQ_WQX-4996680
UTAHDWQ_WQX-4995578
UTAHDWQ_WQX-4995041
WFWQC_UT-4995043

trib_fac_sites <- trib_fac_wq_flow %>% 
  ungroup() %>%
  select(MonitoringLocationIdentifier, MonitoringLocationName, 
           MonitoringLocationTypeName, LatitudeMeasure, LongitudeMeasure) %>%
  distinct() %>%
  mutate(SiteType = ifelse(MonitoringLocationTypeName %in% 
                  c("Canal Drainage", "Canal Irrigation", "Canal Transport"), "Canal", 
                ifelse(MonitoringLocationTypeName %in% c("River/Stream", "Stream", "Stream: Ditch"), "Stream",
                       ifelse(MonitoringLocationTypeName %in% c("Facility Industrial", "Facility Other"), "Facility", 
                              "Storm Sewer"))))


trib_fac_sites_sf <- st_as_sf(trib_fac_sites, coords = c("LongitudeMeasure", "LatitudeMeasure"),
                              crs = 4326, dim = "XY") #CRS = WGS84


trib_fac_sites_intersect <- trib_fac_sites_sf %>%
  st_intersection(wshds) %>%
  left_join(trib_fac_wq_flow, .,) 

unique(trib_fac_sites_intersect$MonitoringLocationName[is.na(trib_fac_sites_intersect$WS_Name)])

##  [1] "Provo River Delta Restoration"               
##  [2] "JORDAN R AT UTAH L OUTLET U121 XING"         
##  [3] "Powell Slough WMA North Outfall to Utah Lake"
##  [4] "Powell Slough WMA South Outfall to Utah Lake"
##  [5] "Dry Ck Near Utah Lake-WLA"                   
##  [6] "Unnamed flow at 4000 West @ mouth"           
##  [7] "American Fork River at mouth"                
##  [8] "Timpanogos WWTP at Utah Lake mouth"          
##  [9] "Lindon Drain at Utah Lake Inlet"             
## [10] "Beer Creek at Utah Lake mouth"               
## [11] "Dry Creek North Fork @ Confluence"           
## [12] "Dry Creek South Fork @ Confluence"           
## [13] "Hobble Creek at Utah Lake mouth"             
## [14] "Mill Race at mouth 1.25 miles west of I-15"

unique(trib_fac_sites_intersect$MonitoringLocationIdentifier[is.na(trib_fac_sites_intersect$WS_Name)])

##  [1] "UTAHDWQ_WQX-4917343" "UTAHDWQ_WQX-4994790" "UTAHDWQ_WQX-4995210"
##  [4] "UTAHDWQ_WQX-4995230" "UTAHDWQ_WQX-4996040" "WFWQC_UT-4917712"   
##  [7] "WFWQC_UT-4994958"    "WFWQC_UT-4995043"    "WFWQC_UT-4995075"   
## [10] "WFWQC_UT-4995210"    "WFWQC_UT-4995467"    "WFWQC_UT-4996040"   
## [13] "WFWQC_UT-4996047"    "WFWQC_UT-4996048"    "WFWQC_UT-4996096"   
## [16] "WFWQC_UT-4996536"

unique(trib_fac_sites_intersect$WS_Name)

##  [1] NA                              "Lehi Drain 2"                 
##  [3] "Tickville Wash"                "Dry Creek - Saratoga"         
##  [5] "Lehi Spring Creek"             "American Fork River"          
##  [7] "Timp SSD"                      "Lindon Drain"                 
##  [9] "Vineyard Drain 6"              "Powell Slough Major"          
## [11] "Currant Creek"                 "Benjamin Slough"              
## [13] "Spanish Fork River"            "Dry Creek - Spanish Fork"     
## [15] "Provo Bay 3"                   "Provo Bay 5"                  
## [17] "Hobble Creek"                  "Spring Creek - Springville"   
## [19] "Mill Race"                     "Provo River"                  
## [21] "4000 South Drain Spanish Fork" "Vineyard Drain 5"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Powell Slough WMA North Outfall to Utah Lake"] <- "Powell Slough Major"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Powell Slough WMA South Outfall to Utah Lake"] <- "Powell Slough Major"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Dry Ck Near Utah Lake-WLA"] <- "Dry Creek - Spanish Fork"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Mill Race at mouth 1.25 miles west of I-15"] <- "Mill Race"  

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Hobble Creek at Utah Lake mouth"] <- "Hobble Creek"
                                 
trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Beer Creek at Utah Lake mouth"] <- "Benjamin Slough"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Lindon Drain at Utah Lake Inlet"] <- "Lindon Drain"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Timpanogos WWTP at Utah Lake mouth"] <- "Timp SSD"
                                   
trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "American Fork River at mouth"] <- "American Fork River"
                                
trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Unnamed flow at 4000 West @ mouth"] <- "4000 South Drain Spanish Fork"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "JORDAN R AT UTAH L OUTLET U121 XING"] <- "Jordan River"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "DRY CK EAST TRIBUTARY"] <- "Dry Creek - Spanish Fork"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "Drainage Canal 0.5 mile bl I-15 at about 2500 West, Springville"] <- "Dry Creek - Spanish Fork"

trib_fac_sites_intersect$WS_Name[trib_fac_sites_intersect$MonitoringLocationName == 
                                   "PROVO STATION 6-WLA"] <- "Mill Race"



trib_fac_sites_wsds <- trib_fac_sites_intersect %>%
  ungroup() %>%
  select(MonitoringLocationIdentifier, OrganizationIdentifier:LongitudeMeasure, 
         SiteType:Shape_Area) %>%
  distinct() %>%
  arrange(WS_Name)
# write.csv(trib_fac_sites_wsds, file = "./Data/Processed/trib_fac_sites_initial.csv",
#           row.names = FALSE)

trib_fac_sites_counts <- trib_fac_sites_intersect %>%
  ungroup() %>%
  group_by(MonitoringLocationIdentifier, OrganizationIdentifier, OrganizationFormalName,
           MonitoringLocationName, MonitoringLocationTypeName, LatitudeMeasure,
           LongitudeMeasure, SiteType, OBJECTID_1, OBJECTID, TARGET_FID,
           Type, wsa_sqkm, WS_Name, Monitored_, 
           Shape_Leng, Shape_Area) %>%
  tally() %>%
  ungroup() %>%
  select(MonitoringLocationIdentifier, n)

# trib_fac_sites_SD <- read.csv("./Data/Processed/trib_fac_sites.csv")
# 
# trib_fac_sites_combined <- full_join(trib_fac_sites_SD, trib_fac_sites_counts)
# write.csv(trib_fac_sites_combined, file = "./Data/Processed/trib_fac_sites_counts.csv",
#           row.names = FALSE)

trib_fac_sites_wsds <- trib_fac_sites_wsds %>%
  full_join(., trib_fac_sites_counts)
# write.csv(trib_fac_sites_wsds, file = "./Data/Processed/trib_fac_sites_counts.csv",
#           row.names = FALSE)

trib_fac_sites_intersect_downstream <- trib_fac_sites_intersect %>%
  filter(MonitoringLocationIdentifier %in% 
           c("UTAHDWQ_WQX-4995312", "UTAHDWQ_WQX-4995310", # Currant Creek, duplicate sites at same lat/lon
             "UTAHDWQ_WQX-4995465", "WFWQC_UT-4995467", # Benjamin Slough
             "UTAHDWQ_WQX-5919910", "WFWQC_UT-4917712", # 4000 South Drain Spanish Fork
             "UTAHDWQ_WQX-4995578", "WFWQC_UT-4995575", # Spanish Fork River
             "UTAHDWQ_WQX-4996040", "WFWQC_UT-4996040", # Dry Creek - Spanish Fork
             "UTAHDWQ_WQX-4996042", "UTAHDWQ_WQX-4996044", # Dry Creek - Spanish Fork
             "WFWQC_UT-4996100", "UTAHDWQ_WQX-4996100", "WFWQC_UT-4996096",
             "UTAHDWQ_WQX-4996275", "WFWQC_UT-4996275", # Hobble Creek
             "WFWQC_UT-4996536", "WFWQC_UT-4996540", "UTAHDWQ_WQX-4996540", "UTAHDWQ_WQX-4996566", # Mill Race
             "WFWQC_UT-4996680", "UTAHDWQ_WQX-4996680", # Provo River
             "UTAHDWQ_WQX-4995230", "WFWQC_UT-4995210", "UTAHDWQ_WQX-4995210", # Powell Slough Major
             "WFWQC_UT-4995075", "UTAHDWQ_WQX-4995120", # Lindon Drain 
             "WFWQC_UT-4995043", "UTAHDWQ_WQX-4995041", "UTAHDWQ_WQX-4995038", # Timp SSD
             "WFWQC_UT-4994958", "UTAHDWQ_WQX-4994960", # American Fork River
             "WFWQC_UT-4994948", "UTAHDWQ_WQX-4994950", #Lehi Spring Creek
             "UTAHDWQ_WQX-4994804", # Dry Creek - Saratoga
             "UTAHDWQ_WQX-4994792")) %>% #Lake Mountain - Israel Canyon
  mutate(Backwater = ifelse(MonitoringLocationIdentifier %in% 
         c("WFWQC_UT-4917712", "WFWQC_UT-4994958", "UTAHDWQ_WQX-4995465", 
           "WFWQC_UT-4995467", "WFWQC_UT-4996536", "UTAHDWQ_WQX-4995210", 
           "UTAHDWQ_WQX-4995230", "WFWQC_UT-4995210", "WFWQC_UT-4995575"), 
         "Yes", "No")) %>%
  left_join(., USGS_flow_processed)
unique(trib_fac_sites_intersect_downstream$WS_Name)

##  [1] "Tickville Wash"                "Dry Creek - Saratoga"         
##  [3] "Lehi Spring Creek"             "American Fork River"          
##  [5] "Timp SSD"                      "Lindon Drain"                 
##  [7] "Powell Slough Major"           "Currant Creek"                
##  [9] "Benjamin Slough"               "Spanish Fork River"           
## [11] "Dry Creek - Spanish Fork"      "Hobble Creek"                 
## [13] "Spring Creek - Springville"    "Mill Race"                    
## [15] "Provo River"                   "4000 South Drain Spanish Fork"

# check that flow at USGS gage is in line with measured flow at chem sites, when data exist together.
# Some cases in Spanish Fork River where measured flow was near zero but USGS flow was not. Otherwise lines up on 1:1 line.
# Will proceed to replace missing flow with USGS flow in those 
# ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name %in% 
#                 c("Hobble Creek", "Spanish Fork River", "Provo River", "Currant Creek")),
#        aes(x = Flow.cfs, y = Flow.USGS.cfs, color = WS_Name)) +
#   geom_point() +
#   geom_abline(slope = 1, intercept = 0)

trib_fac_sites_intersect_downstream$Flow.cfs[is.na(trib_fac_sites_intersect_downstream$Flow.cfs)] <-
  trib_fac_sites_intersect_downstream$Flow.USGS.cfs[is.na(trib_fac_sites_intersect_downstream$Flow.cfs)]

# remove zero flow samples from watersheds that do not go dry


# then, calculate daily loads

trib_fac_sites_intersect_downstream <- trib_fac_sites_intersect_downstream %>%
  mutate(Flow.cfs = case_when(WS_Name == "Lehi Spring Creek" & Flow.cfs == 0 ~ NA_real_, 
                              WS_Name == "Timp SSD" & OrganizationIdentifier == "WFWQC_UT" & 
                                Flow.cfs == 0 ~ NA_real_, 
                              WS_Name == "Powell Slough Major" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Provo River" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Mill Race" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Hobble Creek" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Dry Creek - Spanish Fork" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Spanish Fork River" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "4000 South Drain Spanish Fork" & Flow.cfs == 0 ~ NA_real_,
                              WS_Name == "Benjamin Slough" & Flow.cfs == 0 ~ NA_real_, 
                              TRUE ~ Flow.cfs)) %>%
  mutate(load.TN.kgd = TN * Flow.cfs * 28.3168 * 86400 / 1000000, 
         load.TDN.kgd = TDN * Flow.cfs * 28.3168 * 86400 / 1000000,
         load.TP.kgd = TP * Flow.cfs * 28.3168 * 86400 / 1000000,
         load.TDP.kgd = TDP * Flow.cfs * 28.3168 * 86400 / 1000000,
         load.TOC.kgd = TOC * Flow.cfs * 28.3168 * 86400 / 1000000,
         load.DOC.kgd = DOC * Flow.cfs * 28.3168 * 86400 / 1000000)

trib_fac_sites_intersect_downstream$load.TN.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0
trib_fac_sites_intersect_downstream$load.TDN.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0
trib_fac_sites_intersect_downstream$load.TP.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0
trib_fac_sites_intersect_downstream$load.TDP.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0
trib_fac_sites_intersect_downstream$load.TOC.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0
trib_fac_sites_intersect_downstream$load.DOC.kgd[trib_fac_sites_intersect_downstream$Flow.cfs == 0] <- 0

trib_fac_sites_intersect_downstream$SiteType[is.na(trib_fac_sites_intersect_downstream$SiteType)] <- "Not specified"

trib_fac_sites_downstream <- trib_fac_sites_intersect_downstream %>% 
  ungroup() %>%
  select(MonitoringLocationIdentifier, MonitoringLocationName, OrganizationIdentifier,
           MonitoringLocationTypeName, LatitudeMeasure, LongitudeMeasure, WS_Name) %>%
  distinct()


trib_fac_sites_downstream_sf <- st_as_sf(trib_fac_sites_downstream, coords = c("LongitudeMeasure", "LatitudeMeasure"),
                              crs = 4326, dim = "XY") #CRS = WGS84

trib_fac_sites_downstream_sf$WS_Name <- 
  factor(trib_fac_sites_downstream_sf$WS_Name, 
         levels = c("Tickville Wash", "Dry Creek - Saratoga", "Lehi Spring Creek", 
                    "American Fork River", "Timp SSD", "Lindon Drain", 
                    "Powell Slough Major", "Provo River", "Mill Race", 
                    "Spring Creek - Springville", "Hobble Creek", 
                    "Dry Creek - Spanish Fork", "Spanish Fork River", 
                    "4000 South Drain Spanish Fork", "Benjamin Slough", "Currant Creek"))
                                                                                                                      
trib_fac_sites_intersect_downstream$WS_Name <- 
  factor(trib_fac_sites_intersect_downstream$WS_Name, 
         levels = c("Tickville Wash", "Dry Creek - Saratoga", "Lehi Spring Creek", 
                    "American Fork River", "Timp SSD", "Lindon Drain", 
                    "Powell Slough Major", "Provo River", "Mill Race", 
                    "Spring Creek - Springville", "Hobble Creek", 
                    "Dry Creek - Spanish Fork", "Spanish Fork River", 
                    "4000 South Drain Spanish Fork", "Benjamin Slough", "Currant Creek"))

# Create two datasets: one starting in 2010, one starting in 2015. The later period of record is the main dataset. 

trib_fac_sites_intersect_downstream_early <- trib_fac_sites_intersect_downstream

trib_fac_sites_intersect_downstream <- trib_fac_sites_intersect_downstream %>%
  filter(Year >= 2015)

# convert concentrations to ug/L, not mg/L
# Note: 'Total Kjeldahl nitrogen_Total' was not multiplied by 1000 due to spaces in the name. Variable not being used but will need to be corrected if things change. 
trib_fac_sites_intersect_downstream <- trib_fac_sites_intersect_downstream %>%
  mutate(NH3.d = NH3.d * 1000,
         NH3.t = NH3.t * 1000,
         NO23.d = NO23.d * 1000, 
         NO23.t = NO23.t *1000, 
         TDN = TDN * 1000, 
         TN = TN * 1000, 
         DOC = DOC * 1000,
         TOC = TOC * 1000, 
         TDP = TDP * 1000, 
         TP = TP * 1000, 
         CO3.t = CO3.t * 1000, 
         TKN.t = TKN.t * 1000, 
         TKN.d = TKN.d * 1000, 
         NO3.t = NO3.t * 1000, 
         NO2.t = NO2.t * 1000, 
         PO4 = PO4 * 1000, 
         TKN.t2 = TKN.t2 * 1000)

# There is one TP value from WFWQC of 2 ug/L, which is below the detection limit of 10 ug/L. Set to 1/2 detection limit. 
trib_fac_sites_intersect_downstream$TP[trib_fac_sites_intersect_downstream$TP == 2] <- 5

trib_fac_sites_intersect_downstream_MillRace <- trib_fac_sites_intersect_downstream %>%
  filter(WS_Name == "Mill Race") %>%
  arrange(ActivityStartDate)
# write.csv(trib_fac_sites_intersect_downstream_MillRace, 
#           file = "./Data/Processed/MillRace_WQandFlow.csv", 
#           row.names = FALSE)

# Theron Miller requested flow data. Subset that data to send to him. 
trib_fac_sites_intersect_downstream_flow <- trib_fac_sites_intersect_downstream %>%
  select(MonitoringLocationIdentifier:ResultStatusIdentifier, orig.units:Flow.cfs, WS_Name)
# write.csv(trib_fac_sites_intersect_downstream_flow,
#           file = "./Data/Processed/UL_TributaryFlow_downstreamsites.csv",
#           row.names = FALSE)

# computing proportion of flow from WFWQC, in total and for Lindon Drain and Spanish Fork
trib_fac_sites_intersect_downstream_flow_entities_TN <-
  trib_fac_sites_intersect_downstream %>%
  drop_na(Flow.cfs, TN) %>%
  group_by(OrganizationIdentifier) %>%
  summarise(n = n())
trib_fac_sites_intersect_downstream_flow_entities_TN

## # A tibble: 2 x 2
##   OrganizationIdentifier     n
##   <chr>                  <int>
## 1 UTAHDWQ_WQX              552
## 2 WFWQC_UT                 111

#   OrganizationIdentifier     n
# * <chr>                  <int>
# 1 UTAHDWQ_WQX              552
# 2 WFWQC_UT                 111

trib_fac_sites_intersect_downstream_flow_entities_TP <-
  trib_fac_sites_intersect_downstream %>%
  drop_na(Flow.cfs, TP) %>%
  group_by(OrganizationIdentifier) %>%
  summarise(n = n())
trib_fac_sites_intersect_downstream_flow_entities_TP

## # A tibble: 2 x 2
##   OrganizationIdentifier     n
##   <chr>                  <int>
## 1 UTAHDWQ_WQX              550
## 2 WFWQC_UT                 331

#   OrganizationIdentifier     n
# * <chr>                  <int>
# 1 UTAHDWQ_WQX              550
# 2 WFWQC_UT                 331

trib_fac_sites_intersect_downstream_flow_entities_LDandSF_TP <-
  trib_fac_sites_intersect_downstream %>%
  filter(WS_Name %in% c("Lindon Drain", "Spanish Fork River")) %>%
  drop_na(Flow.cfs, TP) %>%
  group_by(OrganizationIdentifier, WS_Name) %>%
  summarise(n = n())
trib_fac_sites_intersect_downstream_flow_entities_LDandSF_TP

## # A tibble: 4 x 3
## # Groups:   OrganizationIdentifier [2]
##   OrganizationIdentifier WS_Name                n
##   <chr>                  <fct>              <int>
## 1 UTAHDWQ_WQX            Lindon Drain          38
## 2 UTAHDWQ_WQX            Spanish Fork River    31
## 3 WFWQC_UT               Lindon Drain          36
## 4 WFWQC_UT               Spanish Fork River    21

#   OrganizationIdentifier WS_Name                n
#   <chr>                  <fct>              <int>
# 1 UTAHDWQ_WQX            Lindon Drain          38
# 2 UTAHDWQ_WQX            Spanish Fork River    31
# 3 WFWQC_UT               Lindon Drain          36
# 4 WFWQC_UT               Spanish Fork River    21


trib_fac_sites_intersect_downstream_flow_entities_LDandSF_TN <-
  trib_fac_sites_intersect_downstream %>%
  filter(WS_Name %in% c("Lindon Drain", "Spanish Fork River")) %>%
  drop_na(Flow.cfs, TN) %>%
  group_by(OrganizationIdentifier, WS_Name) %>%
  summarise(n = n())
trib_fac_sites_intersect_downstream_flow_entities_LDandSF_TN

## # A tibble: 4 x 3
## # Groups:   OrganizationIdentifier [2]
##   OrganizationIdentifier WS_Name                n
##   <chr>                  <fct>              <int>
## 1 UTAHDWQ_WQX            Lindon Drain          38
## 2 UTAHDWQ_WQX            Spanish Fork River    31
## 3 WFWQC_UT               Lindon Drain          10
## 4 WFWQC_UT               Spanish Fork River     5

#   OrganizationIdentifier WS_Name                n
#   <chr>                  <fct>              <int>
# 1 UTAHDWQ_WQX            Lindon Drain          38
# 2 UTAHDWQ_WQX            Spanish Fork River    31
# 3 WFWQC_UT               Lindon Drain          10
# 4 WFWQC_UT               Spanish Fork River     5


# summarize for report
trib_fac_downstream_summary <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n.2015 = length(ActivityStartDate[Year == 2015]),
            n.2016 = length(ActivityStartDate[Year == 2016]),
            n.2017 = length(ActivityStartDate[Year == 2017]),
            n.2018 = length(ActivityStartDate[Year == 2018]),
            n.2019 = length(ActivityStartDate[Year == 2019]),
            n.2020 = length(ActivityStartDate[Year == 2020]))

trib_fac_downstream_summary_flow <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(Flow.cfs) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(Flow.cfs),
            P10 = quantile(Flow.cfs, 0.10),
            P25 = quantile(Flow.cfs, 0.25),
            Median = quantile(Flow.cfs, 0.5), 
            Mean = mean (Flow.cfs), 
            P75 = quantile(Flow.cfs, 0.75), 
            P90 = quantile(Flow.cfs, 0.90), 
            Max = max(Flow.cfs))

trib_fac_downstream_summary_TP <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(TP) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(TP),
            P10 = quantile(TP, 0.10),
            P25 = quantile(TP, 0.25),
            Median = quantile(TP, 0.5), 
            Mean = mean (TP), 
            P75 = quantile(TP, 0.75), 
            P90 = quantile(TP, 0.90), 
            Max = max(TP))

trib_fac_downstream_summary_TN <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(TN) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(TN),
            P10 = quantile(TN, 0.10),
            P25 = quantile(TN, 0.25),
            Median = quantile(TN, 0.5), 
            Mean = mean (TN), 
            P75 = quantile(TN, 0.75), 
            P90 = quantile(TN, 0.90), 
            Max = max(TN))

trib_fac_downstream_summary_TDP <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(TDP) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(TDP),
            P10 = quantile(TDP, 0.10),
            P25 = quantile(TDP, 0.25),
            Median = quantile(TDP, 0.5), 
            Mean = mean (TDP), 
            P75 = quantile(TDP, 0.75), 
            P90 = quantile(TDP, 0.90), 
            Max = max(TDP))

trib_fac_downstream_summary_TDN <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(TDN) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(TDN),
            P10 = quantile(TDN, 0.10),
            P25 = quantile(TDN, 0.25),
            Median = quantile(TDN, 0.5), 
            Mean = mean (TDN), 
            P75 = quantile(TDN, 0.75), 
            P90 = quantile(TDN, 0.90), 
            Max = max(TDN))

trib_fac_downstream_summary_TOC <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(TOC) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(TOC),
            P10 = quantile(TOC, 0.10),
            P25 = quantile(TOC, 0.25),
            Median = quantile(TOC, 0.5), 
            Mean = mean (TOC), 
            P75 = quantile(TOC, 0.75), 
            P90 = quantile(TOC, 0.90), 
            Max = max(TOC))

trib_fac_downstream_summary_DOC <- trib_fac_sites_intersect_downstream %>%
  group_by(WS_Name, MonitoringLocationIdentifier, OrganizationIdentifier) %>%
  drop_na(DOC) %>%
  summarise(StartDate = min(ActivityStartDate), 
            EndDate = max(ActivityStartDate),
            n = n(), 
            Min = min(DOC),
            P10 = quantile(DOC, 0.10),
            P25 = quantile(DOC, 0.25),
            Median = quantile(DOC, 0.5), 
            Mean = mean (DOC), 
            P75 = quantile(DOC, 0.75), 
            P90 = quantile(DOC, 0.90), 
            Max = max(DOC))

# write.csv(trib_fac_downstream_summary, "./Data/Processed/DownstreamSitesSummary.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_flow, "./Data/Processed/DownstreamSitesSummary_Flow.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_TN, "./Data/Processed/DownstreamSitesSummary_TN.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_TDN, "./Data/Processed/DownstreamSitesSummary_TDN.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_TP, "./Data/Processed/DownstreamSitesSummary_TP.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_TDP, "./Data/Processed/DownstreamSitesSummary_TDP.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_TOC, "./Data/Processed/DownstreamSitesSummary_TOC.csv",
#           row.names = FALSE)
# write.csv(trib_fac_downstream_summary_DOC, "./Data/Processed/DownstreamSitesSummary_DOC.csv",
#           row.names = FALSE)

Jordan River Outflow

jordan_outflow$Date <- as.Date(jordan_outflow$Date, origin = "1899-12-30")
jordan_outflow <- jordan_outflow %>%
  select(Date, Flow)
colnames(jordan_outflow) <- c("Date", "Flow.narrows.cfs")

# if flow at the monitoring site is zero, assign value of zero to narrows. 
jordan <- trib_fac_sites_intersect %>%
  filter(MonitoringLocationIdentifier == "UTAHDWQ_WQX-4994790") %>%
  filter(Year >= 2015) %>%
  left_join(., jordan_outflow, by = c("ActivityStartDate" = "Date")) %>%
  mutate(Flow.narrows.cfs = case_when(Flow.cfs == 0 ~ Flow.cfs,
                                      TRUE ~ Flow.narrows.cfs)) %>%
    mutate(load.TN.kgd = TN * Flow.narrows.cfs * 28.3168 * 86400 / 1000000, 
         load.TDN.kgd = TDN * Flow.narrows.cfs * 28.3168 * 86400 / 1000000,
         load.TP.kgd = TP * Flow.narrows.cfs * 28.3168 * 86400 / 1000000,
         load.TDP.kgd = TDP * Flow.narrows.cfs * 28.3168 * 86400 / 1000000,
         load.TOC.kgd = TOC * Flow.narrows.cfs * 28.3168 * 86400 / 1000000,
         load.DOC.kgd = DOC * Flow.narrows.cfs * 28.3168 * 86400 / 1000000) 

# flow at the narrows is always larger than or equal to the flow measured at the monitoring station. 
ggplot(jordan, aes(x = Flow.cfs, y = Flow.narrows.cfs)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0)

## Warning: Removed 11 rows containing missing values (geom_point).

jordan$load.TN.kgd[jordan$Flow.cfs == 0] <- 0
jordan$load.TDN.kgd[jordan$Flow.cfs == 0] <- 0
jordan$load.TP.kgd[jordan$Flow.cfs == 0] <- 0
jordan$load.TDP.kgd[jordan$Flow.cfs == 0] <- 0
jordan$load.TOC.kgd[jordan$Flow.cfs == 0] <- 0
jordan$load.DOC.kgd[jordan$Flow.cfs == 0] <- 0


loading_jordan <- jordan %>%
  ungroup() %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.narrows.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# temp <- data.frame(approx(loading_jordan$load.TN.kgd, n = 12))
# loading_jordan$load.TN.kgd <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TDN.kgd, n = 12))
# loading_jordan$load.TDN.kgd <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TP.kgd, n = 12))
# loading_jordan$load.TP.kgd <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TDP.kgd, n = 12))
# loading_jordan$load.TDP.kgd <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TOC.kgd, n = 12))
# loading_jordan$load.TOC.kgd <- temp$y
# temp <- data.frame(approx(loading_jordan$load.DOC.kgd, n = 12))
# loading_jordan$load.DOC.kgd <- temp$y
# 
# temp <- data.frame(approx(loading_jordan$load.TN.tonmo, n = 12))
# loading_jordan$load.TN.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TDN.tonmo, n = 12))
# loading_jordan$load.TDN.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TP.tonmo, n = 12))
# loading_jordan$load.TP.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TDP.tonmo, n = 12))
# loading_jordan$load.TDP.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$load.TOC.tonmo, n = 12))
# loading_jordan$load.TOC.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$load.DOC.tonmo, n = 12))
# loading_jordan$load.DOC.tonmo <- temp$y
# temp <- data.frame(approx(loading_jordan$Flow.acftmo, n = 12))
# loading_jordan$Flow.acftmo <- temp$y

loading_jordan_annual <- loading_jordan %>%
  group_by(WS_Name) %>%
  summarise(Flow.cfs = sum(Flow.cfs), 
            load.TN.tonyr = sum(load.TN.tonmo, na.rm = TRUE), 
            load.TDN.tonyr = sum(load.TDN.tonmo, na.rm = TRUE), 
            load.TP.tonyr = sum(load.TP.tonmo, na.rm = TRUE), 
            Flow.acftyr = sum(Flow.acftmo, na.rm = TRUE),
            load.TDP.tonyr = sum(load.TDP.tonmo, na.rm = TRUE), 
            load.TOC.tonyr = sum(load.TOC.tonmo, na.rm = TRUE), 
            load.DOC.tonyr = sum(load.DOC.tonmo, na.rm = TRUE))

# write.csv(loading_jordan, file = "./Data/Processed/loading_Jordan_monthly.csv",
#        row.names = FALSE)
# write.csv(loading_jordan_annual, file = "./Data/Processed/loading_Jordan_annual.csv",
#        row.names = FALSE)

Groundwater (spring and well) data

Precipitation and Evapotranspiration data

DMR data

dmr_updated_concentration <- dmr_updated %>%
  filter(perm_feature_type_desc == "External Outfall") %>%
  filter(monitoring_location_desc == "Effluent Gross") %>%
  filter(parameter_desc %in%
           c("Nitrogen, ammonia total [as N]", "Nitrogen, Kjeldahl, total [as N]",
             "Nitrogen, Kjeldahl, dissolved [as N]", "Nitrite + Nitrate total [as N]",
             "Nitrogen, nitrate total [as N]", "Phosphate, ortho [as P]",
             "Phosphorus, total [as P]", "Flow, in conduit or thru treatment plant")) %>%
  filter(Year >= 2015 & Year <= 2020) %>%
  filter(NPES_Name %in% c("OREM CITY CORP",  
                          "PROVO CITY CORP", "SALEM CITY CORP", "SPANISH FORK CITY CORP", 
                          "SPRINGVILLE- CITY OF", "TIMPANOGOS SPECIAL SERVICE DIS")) %>%
  filter(statistical_base_short_desc %in% c("MX MO AV", "MO MAX", "30DA AVG")) %>%
  select(npdes_id, Year, Month, 
         value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
         NPES_Name, NPDES.Type, NPDESID.Date, parameter_desc) %>%
  # Salem has some flow rates reported as 21 MGD. Their daily max limit is 3 MGD, and all other rates are 0.21. 
  # Changing 21 to 0.21 
  mutate(dmr_value_nmbr = 
           case_when(NPES_Name == "SALEM CITY CORP" & 
                       parameter_desc == "Flow, in conduit or thru treatment plant" ~ 0.21,
                     TRUE ~ dmr_value_nmbr),
         parameter = case_when(parameter_desc == "Nitrogen, ammonia total [as N]" ~ "NH3.mgL",
                               parameter_desc == "Nitrite + Nitrate total [as N]" ~ "NO3.mgL",
                               parameter_desc == "Flow, in conduit or thru treatment plant" ~ "Flow.MGD",
                               parameter_desc == "Phosphorus, total [as P]" ~ "TP.mgL",
                               parameter_desc == "Phosphate, ortho [as P]" ~ "PO4.mgL",
                               parameter_desc == "Nitrogen, Kjeldahl, total [as N]" ~ "TKN.mgL",
                               parameter_desc == "Nitrogen, Kjeldahl, dissolved [as N]" ~ "TKN.mgL")) %>%
  # salem has duplicated month-years that have both a value for nutrients and a zero value. remove the zero value.
  filter(!(NPES_Name == "SALEM CITY CORP" & dmr_value_nmbr == 0))

# PacifiCorps Energy reports flow as DAILY MX or MAXIMUM (always exclusively one or the other)
# Grab this and re-flow with other facilities
dmr_updated_concentration_Lindon <- dmr_updated %>%
  filter(perm_feature_type_desc == "External Outfall") %>%
  filter(monitoring_location_desc == "Effluent Gross") %>%
  filter(parameter_desc %in%
           c("Nitrogen, ammonia total [as N]", "Nitrogen, Kjeldahl, total [as N]",
             "Nitrogen, Kjeldahl, dissolved [as N]", "Nitrite + Nitrate total [as N]",
             "Nitrogen, nitrate total [as N]", "Phosphate, ortho [as P]",
             "Phosphorus, total [as P]", "Flow, in conduit or thru treatment plant")) %>%
  filter(Year >= 2015 & Year <= 2020) %>%
  filter(NPES_Name %in% c("PacifiCorp Energy")) %>%
  #filter(statistical_base_short_desc %in% c("MX MO AV", "MO MAX", "30DA AVG")) %>%
  select(npdes_id, Year, Month, 
         value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
         NPES_Name, NPDES.Type, NPDESID.Date, parameter_desc) %>%
  mutate(parameter = case_when(parameter_desc == "Nitrogen, ammonia total [as N]" ~ "NH3.mgL",
                               parameter_desc == "Nitrite + Nitrate total [as N]" ~ "NO3.mgL",
                               parameter_desc == "Flow, in conduit or thru treatment plant" ~ "Flow.MGD",
                               parameter_desc == "Phosphorus, total [as P]" ~ "TP.mgL",
                               parameter_desc == "Phosphate, ortho [as P]" ~ "PO4.mgL",
                               parameter_desc == "Nitrogen, Kjeldahl, total [as N]" ~ "TKN.mgL",
                               parameter_desc == "Nitrogen, Kjeldahl, dissolved [as N]" ~ "TKN.mgL"))

dmr_updated_concentration <- rbind(dmr_updated_concentration, dmr_updated_concentration_Lindon)

# ggplot(dmr_updated_concentration, 
#        aes(x = parameter_desc, y = dmr_value_nmbr, fill = statistical_base_short_desc)) +
#   geom_boxplot() +
#   theme(legend.position = "top")
# 
# ggplot(subset(dmr_updated_concentration, parameter_desc %in%
#            c("Nitrogen, Kjeldahl, total [as N]",
#              "Nitrogen, Kjeldahl, dissolved [as N]",
#              "Phosphorus, total [as P]", "Flow, in conduit or thru treatment plant",
#              "Flow rate", "Flow, total")), 
#        aes(x = Year, y = dmr_value_nmbr, color = statistical_base_short_desc)) +
#   geom_jitter(alpha = 0.8) +
#   facet_grid(cols = vars(parameter_desc), rows = vars(NPES_Name), scales = "free_y" ) +
#   scale_color_viridis_d(option = "magma", end = 0.8) +
#   theme(legend.position = "top")

dmr_updated_wide <- dmr_updated_concentration %>%
  select(npdes_id, NPES_Name, Year, Month, dmr_value_nmbr, parameter) %>%
  drop_na(dmr_value_nmbr) %>%
  pivot_wider(names_from = parameter, values_from = dmr_value_nmbr) %>%
  mutate(Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.MGD * 31 * 1000000/325851,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.MGD * 30 * 1000000/325851,
                                     Month %in% c(2) ~ Flow.MGD * 28 * 1000000/325851),
         TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ TP.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ TP.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ TP.mgL * Flow.MGD * 28 * 3.78541/1000),
         PO4.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ PO4.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ PO4.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ PO4.mgL * Flow.MGD * 28 * 3.78541/1000),
         TKN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ TKN.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ TKN.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ TKN.mgL * Flow.MGD * 28 * 3.78541/1000),
         NO3.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ NO3.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ NO3.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ NO3.mgL * Flow.MGD * 28 * 3.78541/1000),
         TN.tonmo = TKN.tonmo + NO3.tonmo)
           
# add in Payson data
dmr_updated_Payson <- dmr_Payson %>%
  mutate(npdes_id = NPDES.ID,
         NPES_Name = "PAYSON CITY CORP",
         Year = Year,
         Month = Month,
         Flow.MGD = Effluent.Flow.MGD,
         NH3.mgL = `NH32.mg/L`,
         TP.mgL = `TP.mg/L`,
         PO4.mgL = `Ortho.P.mg/L`,
         TKN.mgL = `TKN.mg/L`,
         NO3.mgL = `Nitrate.mg/L`,
         Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.MGD * 31 * 1000000/325851,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.MGD * 30 * 1000000/325851,
                                     Month %in% c(2) ~ Flow.MGD * 28 * 1000000/325851),
         TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ TP.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ TP.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ TP.mgL * Flow.MGD * 28 * 3.78541/1000),
         PO4.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ PO4.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ PO4.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ PO4.mgL * Flow.MGD * 28 * 3.78541/1000),
         TKN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ TKN.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ TKN.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ TKN.mgL * Flow.MGD * 28 * 3.78541/1000),
         NO3.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ NO3.mgL * Flow.MGD * 31 * 3.78541/1000,
                                     Month %in% c(4, 6, 9, 11) ~ NO3.mgL * Flow.MGD * 30 * 3.78541/1000,
                                     Month %in% c(2) ~ NO3.mgL * Flow.MGD * 28 * 3.78541/1000),
         TN.tonmo = TKN.tonmo + NO3.tonmo) %>%
  select(npdes_id, NPES_Name, Year, Month, Flow.MGD, NH3.mgL, TP.mgL, 
         PO4.mgL, TKN.mgL, NO3.mgL, Flow.acftmo, TP.tonmo, PO4.tonmo, 
         TKN.tonmo, NO3.tonmo, TN.tonmo)

dmr_updated_wide <- rbind(dmr_updated_wide, dmr_updated_Payson)
  
dmr_updated_wide$Flow.acftmo <-  round(dmr_updated_wide$Flow.acftmo, digits = 2)
dmr_updated_wide$TP.tonmo <-  round(dmr_updated_wide$TP.tonmo, digits = 2)
dmr_updated_wide$PO4.tonmo <-  round(dmr_updated_wide$PO4.tonmo, digits = 2)
dmr_updated_wide$TN.tonmo <-  round(dmr_updated_wide$TN.tonmo, digits = 2)


dmr_updated_wide <- dmr_updated_wide %>%
  mutate(MonitoringLocationIdentifier =
           case_when(NPES_Name == "OREM CITY CORP" ~ "UTAHDWQ_WQX-4995250",
                     NPES_Name == "PacifiCorp Energy" ~ "UTAHDWQ_WQX-4995124",
                     NPES_Name == "PAYSON POWER PROJECT" ~ "UTAHDWQ_WQX-4995480",
                     NPES_Name == "PROVO CITY CORP" ~ "UTAHDWQ_WQX-4996560",
                     NPES_Name == "SALEM CITY CORP" ~ "UTAHDWQ_WQX-4995440",
                     NPES_Name == "SPANISH FORK CITY CORP" ~ "UTAHDWQ_WQX-4996020",
                     NPES_Name == "SPRINGVILLE- CITY OF" ~ "UTAHDWQ_WQX-4996280",
                     NPES_Name == "TIMPANOGOS SPECIAL SERVICE DIS" ~ "UTAHDWQ_WQX-4995040"))

dmr_sites <- dmr_updated_wide %>%
  select(NPES_Name, MonitoringLocationIdentifier) %>%
  distinct() %>%
  left_join(., trib_fac_sites) %>%
  filter(!is.na(MonitoringLocationTypeName))

## Joining, by = "MonitoringLocationIdentifier"

dmr_sites_sf <- st_as_sf(dmr_sites, coords = c("LongitudeMeasure", "LatitudeMeasure"),
                              crs = 4326, dim = "XY") #CRS = WGS84
dmr_sites_sf <- st_intersection(dmr_sites_sf, wshds)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#### OLD DATA ####
# dmr_load <- dmr %>%
#   filter(perm_feature_type_desc == "External Outfall") %>%
#   filter(monitoring_location_desc == "Effluent Gross") %>%
#   filter(Nuts.ParmName %in% 
#            c("Nitrogen, ammonia total [as N]", "Nitrogen, Kjeldahl, total [as N]", 
#              "Nitrogen, nitrate total [as N]", "Phosphate, ortho [as P]", 
#              "Phosphorus, total [as P]")) %>%
#   filter(Year >= 2015) %>%
#   select(npdes_id, Year, Month, 
#          value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
#          NPES_Name, NPDES.Type, NPDESID.Date, Nuts.ParmName, Load) %>%
#   na.omit(Load) %>%
#   group_by(npdes_id, Year, Month, 
#            value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
#            NPES_Name, NPDES.Type, NPDESID.Date, Nuts.ParmName) %>%
#   summarise(Load = mean(Load)) %>%
#   ungroup() %>%
#   separate(col = NPDESID.Date, into = c(NA, "Date")) %>%
#   pivot_wider(names_from = "Nuts.ParmName", values_from = "Load")
# 
# dmr_load$Date <- as.Date(as.numeric(dmr_load$Date), origin = "1899-12-30")
# colnames(dmr_load) <- c("NPDES.ID", "Year", "Month", "value.type", "value", 
#                         "units", "stat", "NPDES.Name", "Type", "Date", 
#                         "TP", "NH3", "PO4", "TKN", "NO3")
# 
# dmr_load <- dmr_load %>%
#   filter(stat == "30DA AVG") %>%
#   # load is in ton/yr. convert to metric tons/mo
#   mutate(TP.tonmo = TP/(12*1.10231), 
#          PO4.tonmo = PO4/(12*1.10231),
#          TKN.tonmo = TKN/(12*1.10231),
#          NH3.tonmo = NH3/(12*1.10231),
#          NO3.tonmo = NO3/(12*1.10231)) %>%
#   select(-TP, -PO4, -TKN, -NH3, -NO3) %>%
#   mutate(MonitoringLocationIdentifier = 
#            case_when(NPDES.Name == "MCWANE DUCTILE-UTAH" ~ "UTAHDWQ_WQX-4996430", 
#                      NPDES.Name == "OREM CITY CORP" ~ "UTAHDWQ_WQX-4995250",
#                      NPDES.Name == "PacifiCorp Energy" ~ "UTAHDWQ_WQX-4995124",
#                      NPDES.Name == "PAYSON POWER PROJECT" ~ "UTAHDWQ_WQX-4995480",
#                      NPDES.Name == "PROVO CITY CORP" ~ "UTAHDWQ_WQX-4996560",
#                      NPDES.Name == "SALEM CITY CORP" ~ "UTAHDWQ_WQX-4995440",
#                      NPDES.Name == "SPANISH FORK CITY CORP" ~ "UTAHDWQ_WQX-4996020",
#                      NPDES.Name == "SPRINGVILLE- CITY OF" ~ "UTAHDWQ_WQX-4996280",
#                      NPDES.Name == "TIMPANOGOS SPECIAL SERVICE DIS" ~ "UTAHDWQ_WQX-4995040"))
# 
# 
# dmr_concentration <- dmr %>%
#   filter(perm_feature_type_desc == "External Outfall") %>%
#   filter(Nuts.ParmName %in% 
#            c("Nitrogen, ammonia total [as N]", "Nitrogen, Kjeldahl, total [as N]", 
#              "Nitrogen, nitrate total [as N]", "Phosphate, ortho [as P]", 
#              "Phosphorus, total [as P]")) %>%
#   filter(Year >= 2015) %>%
#   select(npdes_id, Year, Month, 
#          value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
#          NPES_Name, NPDES.Type, NPDESID.Date, Nuts.ParmName, Load) %>%
#   na.omit(Load) %>%
#   group_by(npdes_id, Year, Month, 
#            value_type_desc, dmr_value_nmbr, dmr_unit_desc, statistical_base_short_desc,
#            NPES_Name, NPDES.Type, NPDESID.Date, Nuts.ParmName) %>%
#   summarise(Load = mean(Load)) %>%
#   ungroup() %>%
#   separate(col = NPDESID.Date, into = c(NA, "Date")) %>%
#   filter(dmr_unit_desc == "mg/L" & statistical_base_short_desc == "30DA AVG") %>%
#   pivot_wider(names_from = "Nuts.ParmName", values_from = "dmr_value_nmbr")
# 
# dmr_concentration$Date <- as.Date(as.numeric(dmr_concentration$Date), origin = "1899-12-30")
# colnames(dmr_concentration) <- c("NPDES.ID", "Year", "Month", "value.type", 
#                         "units", "stat", "NPDES.Name", "Type", "Date", "Load",
#                         "PO4", "TP", "NH3",  "TKN")
# 
# dmr_concentration <- dmr_concentration %>%
#   mutate(WS_Name = case_when(NPDES.Name == "TIMPANOGOS SPECIAL SERVICE DIS" ~ "Timp SSD",
#                              NPDES.Name == "OREM CITY CORP" ~ "Powell Slough Major",
#                              NPDES.Name == "PROVO CITY CORP" ~ "Mill Race",
#                              NPDES.Name == "SPRINGVILLE- CITY OF" ~ "Spring Creek - Springville",
#                              NPDES.Name == "SPANISH FORK CITY CORP" ~ "Dry Creek - Spanish Fork", 
#                              NPDES.Name == "SALEM CITY CORP" ~ "Benjamin Slough",
#                              NPDES.Name == "PAYSON CITY CORP" ~ "Benjamin Slough", 
#                              NPDES.Name == "PacifiCorp Energy" ~ "Lindon Drain")) %>%
#   mutate(PO4 = PO4*1000, 
#          TP = TP*1000, 
#          NH3 = NH3*1000, 
#          TKN = TKN*1000)
# 
# 
# dmr_flow <- dmr %>%
#   filter(perm_feature_type_desc == "External Outfall") %>%
#   filter(parameter_desc =="Flow, in conduit or thru treatment plant" & 
#            statistical_base_short_desc == "30DA AVG") %>%
#   filter(Year >= 2015) %>%
#   filter(NPES_Name %in% c("OREM CITY CORP", "PacifiCorp Energy", "PAYSON CITY CORP", 
#                           "PROVO CITY CORP", "SALEM CITY CORP", "SPANISH FORK CITY CORP", 
#                           "SPRINGVILLE- CITY OF", "TIMPANOGOS SPECIAL SERVICE DIS")) %>%
#   select(npdes_id, Year, Month, 
#          dmr_value_standard_units, statistical_base_short_desc,
#          NPES_Name, NPDES.Type, NPDESID.Date) %>%
#   group_by(Month, NPES_Name) %>%
#   summarise(Flow.MGD = mean(dmr_value_standard_units)) 
# 
# # fill in NA values with average of all months for that location
# dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "OREM CITY CORP" & is.na(dmr_flow$Flow.MGD)] <-
#   mean(dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "OREM CITY CORP"], na.rm = TRUE)
# 
# dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "PROVO CITY CORP" & is.na(dmr_flow$Flow.MGD)] <-
#   mean(dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "PROVO CITY CORP"], na.rm = TRUE)
# 
# dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "SALEM CITY CORP" & is.na(dmr_flow$Flow.MGD)] <-
#   mean(dmr_flow$Flow.MGD[dmr_flow$NPES_Name == "SALEM CITY CORP"], na.rm = TRUE)
# 
# dmr_flow <- dmr_flow %>%
#   mutate(Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.MGD * 31 * 3.06889, 
#                                  Month %in% c(4, 6, 9, 11) ~ Flow.MGD * 30 * 3.06889,
#                                  Month %in% c(2) ~ Flow.MGD * 28 * 3.06889)) %>%
#   arrange(NPES_Name, Month)
# 
# dmr_flow_annual <- dmr_flow %>%
#   group_by(NPES_Name) %>%
#   summarise(Flow.acftyr = sum(Flow.acftmo))
# 
# # PacifiCorp Energy does not report 30-day average
# dmr_flow_Lindon <- dmr %>%
#   filter(perm_feature_type_desc == "External Outfall") %>%
#   filter(parameter_desc =="Flow, in conduit or thru treatment plant") %>%
#   filter(Year >= 2015) %>%
#   filter(NPES_Name %in% c("PacifiCorp Energy")) %>%
#   select(npdes_id, Year, Month, 
#          dmr_value_standard_units, statistical_base_short_desc,
#          NPES_Name, NPDES.Type, NPDESID.Date) %>%
#   drop_na() %>%
#   group_by(Year, Month, NPES_Name) %>%
#   summarise(Flow.MGD = round(mean(dmr_value_standard_units), digits = 2)) 
#   
# # # 30 day avg most available for TP
# # ggplot(dmr_load, aes(x = NPDES.Name, y = TP, fill = stat)) +
# #   geom_boxplot()
# # 
# # # 30 day avg most available for PO4
# # ggplot(dmr_load, aes(x = NPDES.Name, y = PO4, fill = stat)) +
# #   geom_boxplot()
# # 
# # # 30 day avg most available for TKN
# # ggplot(dmr_load, aes(x = NPDES.Name, y = TKN, fill = stat)) +
# #   geom_boxplot()
# # 
# # # data very limited for NO3. may want to use TKN instead.
# # ggplot(dmr_load, aes(x = NPDES.Name, y = NO3, fill = stat)) +
# #   geom_boxplot()
# # 
# # # several stats  available for NH3. may want to use TKN instead.
# # ggplot(dmr_load, aes(x = NPDES.Name, y = NH3, fill = stat)) +
# #   geom_boxplot()
# 
# # add in Payson data
# dmr_Payson_load <- dmr_Payson %>%
#   mutate(NPDES.ID = NPDES.ID, 
#          Year = Year,
#          Month = Month, 
#          value.type = NA, 
#          value = NA,
#          units = "mg/L", 
#          stat = "30DA AVG",
#          NPDES.Name = "PAYSON CITY CORP",
#          Type = "Municipal",
#          Date = as.Date(paste(Year, Month, "01", sep = "-")),
#          TP.tonmo = `TP.mg/L`*Effluent.Flow.MGD * 30 * 3785411.78 / 1000000000,
#          PO4.tonmo = `Ortho.P.mg/L`*Effluent.Flow.MGD * 30 * 3785411.78 / 1000000000,
#          TKN.tonmo = `TKN.mg/L`*Effluent.Flow.MGD * 30 * 3785411.78 / 1000000000,
#          NH3.tonmo = NA,
#          NO3.tonmo = `Nitrate.mg/L`*Effluent.Flow.MGD * 30 * 3785411.78 / 1000000000,
#          MonitoringLocationIdentifier = NA) %>%
#   select(-c(Effluent.Flow.MGD, `NH32.mg/L`, `Ortho.P.mg/L`, `TP.mg/L`,
#             `Nitrate.mg/L`, `Nitrite.mg/L`, `TKN.mg/L`))
# 
# dmr_load <- rbind(dmr_load, dmr_Payson_load)
# 
# dmr_Payson_concentration <- dmr_Payson %>%
#   mutate(NPDES.ID = NPDES.ID, 
#          Year = Year,
#          Month = Month, 
#          value.type = NA, 
#          units = "mg/L", 
#          stat = "30DA AVG",
#          NPDES.Name = "PAYSON CITY CORP",
#          Type = "Municipal",
#          Date = as.Date(paste(Year, Month, "01", sep = "-")),
#          Load = NA,
#          PO4 = `Ortho.P.mg/L`*1000,
#          TP = `TP.mg/L`*1000,
#          NH3 = NA,
#          TKN = `TKN.mg/L`*1000,
#          WS_Name = "Benjamin Slough") %>%
#   select(-c(Effluent.Flow.MGD, `NH32.mg/L`, `Ortho.P.mg/L`, `TP.mg/L`,
#             `Nitrate.mg/L`, `Nitrite.mg/L`, `TKN.mg/L`))
# 
# dmr_concentration <- rbind(dmr_concentration, dmr_Payson_concentration)
# 
# # dmr_Payson_flow <- dmr_Payson %>%
# #   mutate(Month = Month, 
# #          NPES_Name = "PAYSON CITY CORP",
# #          Flow.MGD = Effluent.Flow.MGD,
# #          Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.MGD * 31 * 3.06889, 
# #                                  Month %in% c(4, 6, 9, 11) ~ Flow.MGD * 30 * 3.06889,
# #                                  Month %in% c(2) ~ Flow.MGD * 28 * 3.06889)) %>%
# #   select(-c(NPDES.ID, Year, Effluent.Flow.MGD, `NH32.mg/L`, `Ortho.P.mg/L`,
# #             `TP.mg/L`, `Nitrate.mg/L`, `Nitrite.mg/L`, `TKN.mg/L`))
# #   
# # dmr_flow <- rbind(dmr_flow, dmr_Payson_flow)
# 
# dmr_sites <- dmr_load %>%
#   select(NPDES.Name, MonitoringLocationIdentifier) %>%
#   distinct() %>%
#   left_join(., trib_fac_sites) %>%
#   filter(!is.na(MonitoringLocationTypeName))
# 
# dmr_sites_sf <- st_as_sf(dmr_sites, coords = c("LongitudeMeasure", "LatitudeMeasure"),
#                               crs = 4326, dim = "XY") #CRS = WGS84
# dmr_sites_sf <- st_intersection(dmr_sites_sf, wshds)
# 
# 
# 
# #write.csv(dmr_flow, file = "./Data/Processed/DMR_flow.csv", row.names = FALSE)

Visualization

Maps

The top graph depicts all tributary and facility monitoring sites. The second graph depicts each monitored watershed individually from the NW side of the lake moving clockwise, with downstream monitoring sites marked. The third graph depicts all watersheds with only downstream monitoring sites marked.

Water Quality Time Series

Flow time series

Load time series

Note: this is an exploration, not a quantification of total load. Several watersheds have >2 sites that make up total load.

# examine Mill race, specifically
ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Mill Race"),
       aes(x = ActivityStartDate, y = load.TN.kgd, 
           color = OrganizationIdentifier)) +
  geom_point(alpha = 0.7, size = 2) +
  scale_y_log10() +
  labs(x = "", y = "TN load (kg/d)", shape = "Site Type", color = "Organization") +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(ncol = 2))

## Warning: Removed 72 rows containing missing values (geom_point).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Mill Race"),
       aes(x = ActivityStartDate, y = load.TP.kgd, 
           color = OrganizationIdentifier)) +
  geom_point(alpha = 0.7, size = 2) +
  scale_y_log10() +
  labs(x = "", y = "TP load (kg/d)", shape = "Site Type", color = "Organization") +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(ncol = 2))

## Warning: Removed 41 rows containing missing values (geom_point).

# ggplot(dmr_load, aes(x = Date, y = TP.tonmo)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TP load (metric ton/mo)")
# 
# ggplot(dmr_load, aes(x = Date, y = PO4.tonmo)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TDP load (metric ton/mo)")
# 
# ggplot(dmr_load, aes(x = Date, y = TKN.tonmo)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TKN load (metric ton/mo)")
# 
# 
# ggplot(dmr_load, aes(x = as.factor(Month), y = TP.tonmo)) +
#   geom_boxplot() +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TP load (metric ton/mo)", x = "Month")
# 
# ggplot(dmr_load, aes(x = as.factor(Month), y = PO4.tonmo)) +
#   geom_boxplot() +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TDP load (metric ton/mo)", x = "Month")
# 
# ggplot(dmr_load, aes(x = as.factor(Month), y = TKN.tonmo)) +
#   geom_boxplot() +
#   facet_wrap(vars(NPDES.Name), nrow = 2) +
#   scale_y_log10() +
#   labs(y = "TKN load (metric ton/mo)", x = "Month")

DMR Flow, concentration, and load

ggplot(dmr_updated_wide, aes(x = Month, y = Flow.acftmo)) +
  geom_line() +
  labs(y = "Flow (ac*ft/mo)") +
  scale_y_continuous(sec.axis = sec_axis(~./(30 * 3.06889), name = "Flow (MGD)")) +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

#ggsave("./Code/Output/DMR_Flow_Monthly.jpg", width = 8, height = 4)

ggplot(dmr_updated_wide, aes(x = lubridate::mdy(paste(Month, 1, Year)), y = TP.mgL)) +
  geom_point() +
  scale_y_log10() +
  labs(y = expression("TP (mg/L)"), x = "") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 78 rows containing missing values (geom_point).

#ggsave("./Code/Output/DMR_TPconc_TimeSeries.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = lubridate::mdy(paste(Month, 1, Year)), y = TKN.mgL)) +
  geom_point() +
  scale_y_log10() +
  labs(y = expression("TKN (mg/L)"), x = "") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 188 rows containing missing values (geom_point).

#ggsave("./Code/Output/DMR_TKNconc_TimeSeries.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = lubridate::mdy(paste(Month, 1, Year)), y = NO3.mgL)) +
  geom_point() +
  scale_y_log10() +
  labs(y = expression("Nitrate (mg/L)"), x = "") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 209 rows containing missing values (geom_point).

#ggsave("./Code/Output/DMR_NO3conc_TimeSeries.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = lubridate::mdy(paste(Month, 1, Year)), y = PO4.mgL)) +
  geom_point() +
  scale_y_log10() +
  labs(y = expression("Orthophosphate (mg/L)"), x = "") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 200 rows containing missing values (geom_point).

#ggsave("./Code/Output/DMR_PO4conc_TimeSeries.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = TP.mgL)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = expression("TP (mg/L)"), x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 78 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TPconc_Monthly.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Year), y = TP.mgL)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = expression("TP (mg/L)"), x = "Year") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 78 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TPconc_annual.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = TKN.mgL)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = expression("TKN (mg/L)"), x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 188 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TKNconc_Monthly.jpg", width = 8, height = 5)


ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = NO3.mgL)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = expression("Nitrate (mg/L)"), x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 209 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TKNconc_Monthly.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = PO4.mgL)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = expression("Orthophosphate (mg/L)"), x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 201 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_PO4conc_Monthly.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = TP.tonmo)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = "TP (metric ton/mo)", x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 81 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TPload_Monthly.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = TN.tonmo)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = "TN load (metric ton/mo)", x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Removed 210 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_TNload_Monthly.jpg", width = 8, height = 5)

ggplot(dmr_updated_wide, aes(x = as.factor(Month), y = PO4.tonmo)) +
  geom_boxplot() +
  scale_y_log10() +
  labs(y = "Orthophosphate load (metric ton/mo)", x = "Month") +
  facet_wrap(vars(NPES_Name), nrow = 2,  
             labeller = label_wrap_gen(width = 20, multi_line = TRUE))

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 203 rows containing non-finite values (stat_boxplot).

#ggsave("./Code/Output/DMR_PO4load_Monthly.jpg", width = 8, height = 5)

Monthly loads by watershed

Tickville Wash

There is only one inflow monitored in Tickville Wash. Estimates are based on this site.

Dry Creek - Saratoga

There is only one inflow monitored in Dry Creek - Saratoga. Estimates are based on this site.

Lehi Spring Creek

There are two sites in Lehi Spring Creek. DWQ site is upstream of WFWQC site. Will compare loads to determine if the estimates can be combined.

American Fork River

There are two sites in American Fork River. DWQ site is upstream of WFWQC site. Will compare loads to determine if the estimates can be combined.

Timp SSD

REVISIT East side tributary (UTAHDWQ_WQX-4995041) has measurable concentrations but no flow (one near-zero flow, one flow duplicated from 4995038). Scott’s recommendation was to average concentrations at 4995041 and 4995038 and apply them to the flow at 4995038, but this assumes that the volume-weighted average is the same as the arithmetic average (i.e., that the volumes at both sites are the same). I’m not sure if this assumption is valid. For now, I will ignore the east side tributary and calculate loads based on the sites that have flow data.

Lindon Drain

There are two sites in Lindon Drain. DWQ site is upstream of WFWQC site. DWQ site needs to be added to DMR data to represent WFWQC loading (WFWQC site represents total loading, DWQ site is located upstream of power plant effluent).

Powell Slough Major

Powell Slough has both North and South outfalls. Need to add together load from North and South to compute total load. WFWQC samples only at the north outfall.

Provo River

There are two site codes in Provo River (one DWQ and one WFWQC), but these are in the same location. Will compare loads to determine if the estimates can be combined.

Mill Race

4996536 represents the sum total of loading from the watershed. 4996540 and 4996566 should be added together to represent total flow.

Spring Creek - Springville

4996275 represents Springville WWTP.

Hobble Creek

4996096 and 4996100 are equivalent sites.

Dry Creek - Spanish Fork

4996040 is the primary site in this watershed and represents the sum total of loads. 4996042 and 4996044 represent equivalent sites for the East tributary but do not include the south fork. East tributary sites will be removed and the primary site used for loads.

Spanish Fork River

There are two sites in Spanish Fork River. DWQ site is upstream of WFWQC site. Will compare loads to determine if the estimates can be combined.

4000 South Drain Spanish Fork

There are two sites in 4000 South Drain Spanish Fork. DWQ site is upstream of WFWQC site. Will compare loads to determine if the estimates can be combined.

Benjamin Slough

There are two sites in Benjamin Slough. DWQ site is upstream of WFWQC site. Will compare loads to determine if the estimates can be combined.

Currant Creek

Two sites, UTAHDWQ_WQX-4995310 and UTAHDWQ_WQX-4995312, are located at the same coordinates. These sites are pooled in load calculations, as they represent the same location.

Monthly and annual load across watersheds

Monitored tributaries

# Tickville Wash
loading_TickvilleWash <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Tickville Wash") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Dry Creek - Saratoga
loading_DryCreekSaratoga <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Dry Creek - Saratoga") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Lehi Spring Creek
# combine data from DWQ and WFWQC data from two sites
# TP and TDP concentrations are below reporting limit for WFWQC --> use only DWQ samples
loading_LehiSpringCreek <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Lehi Spring Creek") %>%
  mutate(load.TP.kgd = case_when(OrganizationIdentifier == "WFWQC_UT" ~ NA_real_,
                                 TRUE ~ load.TP.kgd),
         load.TDP.kgd = case_when(OrganizationIdentifier == "WFWQC_UT" ~ NA_real_,
                                 TRUE ~ load.TDP.kgd)) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# American Fork River
# combine data from DWQ and WFWQC data from two sites
# P and N concentrations are below reporting limit for WFWQC --> use only DWQ samples

loading_AmericanForkRiver <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "American Fork River") %>%
  filter(OrganizationIdentifier != "WFWQC_UT") %>%
  group_by(WS_Name,  Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Timp SSD
# combine data from DWQ and WFWQC data from two sites
loading_TimpSSD <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Timp SSD") %>%
  filter(MonitoringLocationIdentifier != "UTAHDWQ_WQX-4995041") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Lindon Drain
# WFWQC represents most downstream flow, DWQ data should be added to the DMR data for PacifiCorps Energy to represent total flow and load.

loading_LindonDrain_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "PacifiCorp Energy") %>%
  group_by(NPES_Name, Month) %>%
  summarise(load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE),
            Flow.MGD = mean(Flow.MGD, na.rm = TRUE),
            Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE)) %>%
  distinct() %>%
  ungroup() %>%
  mutate(WS_Name = "Lindon Drain", 
         Month = Month,
         Flow.cfs.dmr = Flow.MGD/0.646317,
         Flow.acftmo.dmr = Flow.acftmo,
         load.TP.tonmo.dmr = load.TP.tonmo) %>%
  select(-Flow.MGD, -load.TP.tonmo)

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

loading_LindonDrain_DWQ <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Lindon Drain") %>%
  filter(OrganizationIdentifier == "UTAHDWQ_WQX") %>%
  group_by(WS_Name, Month, ActivityStartDate) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct() %>%
# add average monthly DMR values to DWQ values
  left_join(., loading_LindonDrain_dmr) %>%
  mutate(Flow.cfs = Flow.cfs + Flow.cfs.dmr, 
         Flow.acftmo = Flow.acftmo + Flow.acftmo.dmr, 
         load.TP.tonmo = load.TP.tonmo + load.TP.tonmo.dmr) %>%
  select(-NPES_Name, -Flow.cfs.dmr, -Flow.acftmo.dmr, -load.TP.tonmo.dmr)

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

## Joining, by = c("WS_Name", "Month", "Flow.acftmo")

loading_LindonDrain_WFWQC <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Lindon Drain") %>%
  filter(OrganizationIdentifier == "WFWQC_UT") %>%
  group_by(WS_Name, ActivityStartDate, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE),
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'ActivityStartDate'. You can override using the `.groups` argument.

# combine DWQ and WFWQC data, calculate average 
loading_LindonDrain <- loading_LindonDrain_DWQ %>%
  rbind(., loading_LindonDrain_WFWQC) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE),
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TN.tonmo = mean(load.TN.tonmo, na.rm = TRUE),
            load.TDN.tonmo = mean(load.TDN.tonmo, na.rm = TRUE),
            load.TP.tonmo = mean(load.TP.tonmo, na.rm = TRUE),
            load.TDP.tonmo = mean(load.TDP.tonmo, na.rm = TRUE),
            load.TOC.tonmo = mean(load.TOC.tonmo, na.rm = TRUE),
            load.DOC.tonmo = mean(load.DOC.tonmo, na.rm = TRUE))

## `summarise()` has grouped output by 'WS_Name'. You can override using the `.groups` argument.

# Powell Slough 
# North Outfall
loading_PowellSloughMajor_north <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Powell Slough Major") %>%
  filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4995210", "WFWQC_UT-4995210")) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# interpolate missing months

temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TN.kgd, n = 12))
loading_PowellSloughMajor_north$load.TN.kgd <- temp$y
temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TP.kgd, n = 12))
loading_PowellSloughMajor_north$load.TP.kgd <- temp$y
temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TDP.kgd, n = 12))
loading_PowellSloughMajor_north$load.TDP.kgd <- temp$y
temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TN.tonmo, n = 12))
loading_PowellSloughMajor_north$load.TN.tonmo <- temp$y
temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TP.tonmo, n = 12))
loading_PowellSloughMajor_north$load.TP.tonmo <- temp$y
temp <- data.frame(approx(loading_PowellSloughMajor_north$load.TDP.tonmo, n = 12))
loading_PowellSloughMajor_north$load.TDP.tonmo <- temp$y

loading_PowellSloughMajor_north$load.TDN.kgd[is.na(loading_PowellSloughMajor_north$load.TDN.kgd)] <-
  mean(c(loading_PowellSloughMajor_north$load.TDN.kgd[4], loading_PowellSloughMajor_north$load.TDN.kgd[12]))
loading_PowellSloughMajor_north$load.TOC.kgd[is.na(loading_PowellSloughMajor_north$load.TOC.kgd)] <-
  mean(c(loading_PowellSloughMajor_north$load.TOC.kgd[4], loading_PowellSloughMajor_north$load.TOC.kgd[12]))
loading_PowellSloughMajor_north$load.DOC.tonmo[is.na(loading_PowellSloughMajor_north$load.DOC.tonmo)] <-
  mean(c(loading_PowellSloughMajor_north$load.DOC.tonmo[4], loading_PowellSloughMajor_north$load.DOC.tonmo[12]))
loading_PowellSloughMajor_north$load.TDN.tonmo[is.na(loading_PowellSloughMajor_north$load.TDN.tonmo)] <-
  mean(c(loading_PowellSloughMajor_north$load.TDN.tonmo[4], loading_PowellSloughMajor_north$load.TDN.tonmo[12]))
loading_PowellSloughMajor_north$load.TOC.tonmo[is.na(loading_PowellSloughMajor_north$load.TOC.tonmo)] <-
  mean(c(loading_PowellSloughMajor_north$load.TOC.tonmo[4], loading_PowellSloughMajor_north$load.TOC.tonmo[12]))
loading_PowellSloughMajor_north$load.DOC.tonmo[is.na(loading_PowellSloughMajor_north$load.DOC.tonmo)] <-
  mean(c(loading_PowellSloughMajor_north$load.DOC.tonmo[4], loading_PowellSloughMajor_north$load.DOC.tonmo[12]))


# South Outfall
loading_PowellSloughMajor_south <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Powell Slough Major") %>%
  filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4995230")) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Add together north and south outfalls
loading_PowellSloughMajor <- loading_PowellSloughMajor_north %>%
  rbind(., loading_PowellSloughMajor_south) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = sum(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = sum(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = sum(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = sum(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = sum(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = sum(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = sum(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = sum(Flow.acftmo, na.rm = TRUE),
            load.TN.tonmo = sum(load.TN.tonmo, na.rm = TRUE),
            load.TDN.tonmo = sum(load.TDN.tonmo, na.rm = TRUE),
            load.TP.tonmo = sum(load.TP.tonmo, na.rm = TRUE),
            load.TDP.tonmo = sum(load.TDP.tonmo, na.rm = TRUE),
            load.TOC.tonmo = sum(load.TOC.tonmo, na.rm = TRUE),
            load.DOC.tonmo = sum(load.DOC.tonmo, na.rm = TRUE))

## `summarise()` has grouped output by 'WS_Name'. You can override using the `.groups` argument.

# loading_PowellSloughMajor_initial <- trib_fac_sites_intersect_downstream %>%
#   ungroup() %>%
#   filter(WS_Name == "Powell Slough Major") %>%
#   group_by(WS_Name, OrganizationIdentifier, MonitoringLocationIdentifier, Month) %>%
#   summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
#             load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
#             load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
#             load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
#             load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
#             load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
#             load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
#             load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.TN.kgd * 28/1000),
#             load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.TDN.kgd * 28/1000),
#             load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.TP.kgd * 28/1000),
#             load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.TDP.kgd * 28/1000),
#             load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.TOC.kgd * 28/1000),
#             load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
#                                      Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
#                                      Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
#   distinct() 

# calculate load from N and S outfalls with WFWQC data from N outfall (DWQ site only at S outfall)
# loading_PowellSloughMajor_WFWQCloads <- loading_PowellSloughMajor_initial %>%
#   filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4995230", "WFWQC_UT-4995210")) %>%
#   group_by(WS_Name, Month) %>%
#   summarise(Flow.cfs = sum(Flow.cfs, na.rm = TRUE),
#             load.TN.kgd = sum(load.TN.kgd, na.rm = TRUE), 
#             load.TDN.kgd = sum(load.TDN.kgd, na.rm = TRUE),
#             load.TP.kgd = sum(load.TP.kgd, na.rm = TRUE),
#             load.TDP.kgd = sum(load.TDP.kgd, na.rm = TRUE),
#             load.TOC.kgd = sum(load.TOC.kgd, na.rm = TRUE),
#             load.DOC.kgd = sum(load.DOC.kgd, na.rm = TRUE),
#             load.TN.tonmo = sum(load.TN.tonmo, na.rm = TRUE),
#             load.TDN.tonmo = sum(load.TDN.tonmo, na.rm = TRUE),
#             load.TP.tonmo = sum(load.TP.tonmo, na.rm = TRUE),
#             load.TDP.tonmo = sum(load.TDP.tonmo, na.rm = TRUE),
#             load.TOC.tonmo = sum(load.TOC.tonmo, na.rm = TRUE),
#             load.DOC.tonmo = sum(load.DOC.tonmo, na.rm = TRUE)) %>%
#   mutate(MonitoringLocationIdentifier = "WFWQC_UT-4995210 + UTAHDWQ_WQX-4995230", 
#          OrganizationIdentifier = "WFWQC_UT") %>%
#   select(WS_Name, OrganizationIdentifier, MonitoringLocationIdentifier, Month, 
#          Flow.cfs, load.TN.kgd, load.TDN.kgd, load.TP.kgd, load.TDP.kgd, 
#          load.TOC.kgd, load.DOC.kgd, load.TN.tonmo, load.TDN.tonmo, load.TP.tonmo,
#          load.TDP.tonmo, load.TOC.tonmo, load.DOC.tonmo)

# calculate load from N and S outfalls with DWQ data 
# loading_PowellSloughMajor_UTAHDWQloads <- loading_PowellSloughMajor_initial %>%
#   filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4995230", "UTAHDWQ_WQX-4995210")) %>%
#   group_by(WS_Name, Month) %>%
#   summarise(Flow.cfs = sum(Flow.cfs, na.rm = TRUE),
#             load.TN.kgd = sum(load.TN.kgd, na.rm = TRUE), 
#             load.TDN.kgd = sum(load.TDN.kgd, na.rm = TRUE),
#             load.TP.kgd = sum(load.TP.kgd, na.rm = TRUE),
#             load.TDP.kgd = sum(load.TDP.kgd, na.rm = TRUE),
#             load.TOC.kgd = sum(load.TOC.kgd, na.rm = TRUE),
#             load.DOC.kgd = sum(load.DOC.kgd, na.rm = TRUE),
#             load.TN.tonmo = sum(load.TN.tonmo, na.rm = TRUE),
#             load.TDN.tonmo = sum(load.TDN.tonmo, na.rm = TRUE),
#             load.TP.tonmo = sum(load.TP.tonmo, na.rm = TRUE),
#             load.TDP.tonmo = sum(load.TDP.tonmo, na.rm = TRUE),
#             load.TOC.tonmo = sum(load.TOC.tonmo, na.rm = TRUE),
#             load.DOC.tonmo = sum(load.DOC.tonmo, na.rm = TRUE)) %>%
#   mutate(MonitoringLocationIdentifier = "UTAHDWQ_WQX-4995210 + UTAHDWQ_WQX-4995230", 
#          OrganizationIdentifier = "UTAHDWQ_WQX") %>%
#   select(WS_Name, OrganizationIdentifier, MonitoringLocationIdentifier, Month, 
#          Flow.cfs, load.TN.kgd, load.TDN.kgd, load.TP.kgd, load.TDP.kgd, 
#          load.TOC.kgd, load.DOC.kgd, load.TN.tonmo, load.TDN.tonmo, load.TP.tonmo,
#          load.TDP.tonmo, load.TOC.tonmo, load.DOC.tonmo)

#combine loading from different organizations
#loading_PowellSloughMajor <- rbind(loading_PowellSloughMajor_UTAHDWQloads, loading_PowellSloughMajor_WFWQCloads)

# Provo River
# P and N concentrations are below reporting limit for WFWQC --> use only DWQ samples

loading_ProvoRiver <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Provo River") %>%
  filter(OrganizationIdentifier != "WFWQC_UT") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Mill Race
loading_MillRace_initial <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Mill Race") %>%
  # leave out WFWQC_UT-4996536. Only 3 samples available.
  filter(MonitoringLocationIdentifier != "WFWQC_UT-4996536") %>%
  group_by(WS_Name, OrganizationIdentifier, MonitoringLocationIdentifier, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'OrganizationIdentifier', 'MonitoringLocationIdentifier', 'Month'. You can override using the `.groups` argument.

# separate 4996540. This site represents a portion of the total load, which will be added to 4996566.

loading_MillRace_4996540 <- loading_MillRace_initial %>%
  filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4996540", "WFWQC_UT-4996540")) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TN.tonmo = mean(load.TN.tonmo, na.rm = TRUE),
            load.TDN.tonmo = mean(load.TDN.tonmo, na.rm = TRUE),
            load.TP.tonmo = mean(load.TP.tonmo, na.rm = TRUE),
            load.TDP.tonmo = mean(load.TDP.tonmo, na.rm = TRUE),
            load.TOC.tonmo = mean(load.TOC.tonmo, na.rm = TRUE),
            load.DOC.tonmo = mean(load.DOC.tonmo, na.rm = TRUE)) %>%
  mutate(MonitoringLocationIdentifier = "4996540", 
         OrganizationIdentifier = "UTAHDWQ_WQX") %>%
  select(WS_Name, OrganizationIdentifier, MonitoringLocationIdentifier, Month, 
         Flow.cfs, load.TN.kgd, load.TDN.kgd, load.TP.kgd, load.TDP.kgd, 
         load.TOC.kgd, load.DOC.kgd, Flow.acftmo, load.TN.tonmo, load.TDN.tonmo, 
         load.TP.tonmo, load.TDP.tonmo, load.TOC.tonmo, load.DOC.tonmo)

## `summarise()` has grouped output by 'WS_Name'. You can override using the `.groups` argument.

# calculate load from 4996566
loading_MillRace_4996566 <- loading_MillRace_initial %>%
  filter(MonitoringLocationIdentifier %in% c("UTAHDWQ_WQX-4996566"))

# calculate sum of 4996540 and 4996566 These represent the total load.
loading_MillRace <- loading_MillRace_4996540 %>%
  rbind(loading_MillRace_4996566) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = sum(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = sum(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = sum(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = sum(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = sum(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = sum(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = sum(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = sum(Flow.acftmo, na.rm = TRUE),
            load.TN.tonmo = sum(load.TN.tonmo, na.rm = TRUE),
            load.TDN.tonmo = sum(load.TDN.tonmo, na.rm = TRUE),
            load.TP.tonmo = sum(load.TP.tonmo, na.rm = TRUE),
            load.TDP.tonmo = sum(load.TDP.tonmo, na.rm = TRUE),
            load.TOC.tonmo = sum(load.TOC.tonmo, na.rm = TRUE),
            load.DOC.tonmo = sum(load.DOC.tonmo, na.rm = TRUE))

## `summarise()` has grouped output by 'WS_Name'. You can override using the `.groups` argument.

# Spring Creek - Springville
# combine data from DWQ and WFWQC data from two sites
loading_SpringCreek <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Spring Creek - Springville") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# test differences in DWQ and WFWQC
loading_SpringCreek_DWQ <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(OrganizationIdentifier != "WFWQC_UT") %>%
  filter(WS_Name == "Spring Creek - Springville") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

loading_SpringCreek_WFWQC <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(OrganizationIdentifier == "WFWQC_UT") %>%
  filter(WS_Name == "Spring Creek - Springville") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Hobble Creek
# combine data from DWQ and WFWQC data from three sites - all are equivalent
# TP and TDP concentrations are below reporting limit for WFWQC --> use only DWQ samples
 
loading_HobbleCreek <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Hobble Creek") %>%
  mutate(load.TP.kgd = case_when(OrganizationIdentifier == "WFWQC_UT" ~ NA_real_,
                                 TRUE ~ load.TP.kgd),
         load.TDP.kgd = case_when(OrganizationIdentifier == "WFWQC_UT" ~ NA_real_,
                                 TRUE ~ load.TDP.kgd)) %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Dry Creek - Spanish Fork
loading_DryCreekSpanishFork <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Dry Creek - Spanish Fork") %>%
  filter(MonitoringLocationIdentifier != "UTAHDWQ_WQX-4996042" &
           MonitoringLocationIdentifier != "UTAHDWQ_WQX-4996044") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

loading_DryCreekSpanishFork_January <- 
  data.frame(WS_Name = "Dry Creek - Spanish Fork", Month = 1, 
             Flow.cfs = mean(c(loading_DryCreekSpanishFork$Flow.cfs[1], loading_DryCreekSpanishFork$Flow.cfs[11])),
             load.TN.kgd = mean(c(loading_DryCreekSpanishFork$load.TN.kgd[1], loading_DryCreekSpanishFork$load.TN.kgd[11])),
             load.TDN.kgd = mean(c(loading_DryCreekSpanishFork$load.TDN.kgd[1], loading_DryCreekSpanishFork$load.TDN.kgd[11])),
             load.TP.kgd = mean(c(loading_DryCreekSpanishFork$load.TP.kgd[1], loading_DryCreekSpanishFork$load.TP.kgd[11])),
             load.TDP.kgd = mean(c(loading_DryCreekSpanishFork$load.TDP.kgd[2], loading_DryCreekSpanishFork$load.TDP.kgd[11])),
             load.TOC.kgd = mean(c(loading_DryCreekSpanishFork$load.TOC.kgd[1], loading_DryCreekSpanishFork$load.TOC.kgd[11])),
             load.DOC.kgd = mean(c(loading_DryCreekSpanishFork$load.DOC.kgd[1], loading_DryCreekSpanishFork$load.DOC.kgd[11])),
             Flow.acftmo = mean(c(loading_DryCreekSpanishFork$Flow.acftmo[1], loading_DryCreekSpanishFork$Flow.acftmo[11])),
             load.TN.tonmo = mean(c(loading_DryCreekSpanishFork$load.TN.tonmo[1], loading_DryCreekSpanishFork$load.TN.tonmo[11])),
             load.TDN.tonmo = mean(c(loading_DryCreekSpanishFork$load.TDN.tonmo[1], loading_DryCreekSpanishFork$load.TDN.tonmo[11])),
             load.TP.tonmo = mean(c(loading_DryCreekSpanishFork$load.TP.tonmo[1], loading_DryCreekSpanishFork$load.TP.tonmo[11])),
             load.TDP.tonmo = mean(c(loading_DryCreekSpanishFork$load.TDP.tonmo[2], loading_DryCreekSpanishFork$load.TDP.tonmo[11])),
             load.TOC.tonmo = mean(c(loading_DryCreekSpanishFork$load.TOC.tonmo[1], loading_DryCreekSpanishFork$load.TOC.tonmo[11])),
             load.DOC.tonmo = mean(c(loading_DryCreekSpanishFork$load.DOC.tonmo[1], loading_DryCreekSpanishFork$load.DOC.tonmo[11])))

loading_DryCreekSpanishFork <- rbind(loading_DryCreekSpanishFork_January, loading_DryCreekSpanishFork)

loading_DryCreekSpanishFork$load.TDP.kgd[is.na(loading_DryCreekSpanishFork$load.TDP.kgd)] <-
  mean(c(loading_DryCreekSpanishFork$load.TDP.kgd[3], loading_DryCreekSpanishFork$load.TDP.kgd[12]))
loading_DryCreekSpanishFork$load.TDP.tonmo[is.na(loading_DryCreekSpanishFork$load.TDP.tonmo)] <-
  mean(c(loading_DryCreekSpanishFork$load.TDP.tonmo[3], loading_DryCreekSpanishFork$load.TDP.tonmo[12]))

# Spanish Fork

loading_SpanishFork <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Spanish Fork River") %>%
  group_by(WS_Name,  Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# 4000 South Drain Spanish Fork
# combine data from DWQ and WFWQC data from two sites

loading_4000SouthDrain <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "4000 South Drain Spanish Fork") %>%
  group_by(WS_Name,  Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

#interpolate for missing months
temp <- data.frame(approx(loading_4000SouthDrain$load.TN.kgd, n = 12))
loading_4000SouthDrain$load.TN.kgd <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.TDN.kgd, n = 12))
loading_4000SouthDrain$load.TDN.kgd <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.TOC.kgd, n = 12))
loading_4000SouthDrain$load.TOC.kgd <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.DOC.kgd, n = 12))
loading_4000SouthDrain$load.DOC.kgd <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.TN.tonmo, n = 12))
loading_4000SouthDrain$load.TN.tonmo <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.TDN.tonmo, n = 12))
loading_4000SouthDrain$load.TDN.tonmo <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.TOC.tonmo, n = 12))
loading_4000SouthDrain$load.TOC.tonmo <- temp$y
temp <- data.frame(approx(loading_4000SouthDrain$load.DOC.tonmo, n = 12))
loading_4000SouthDrain$load.DOC.tonmo <- temp$y


# Benjamin Slough
# combine data from DWQ and WFWQC data from two sites

loading_BenjaminSlough <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Benjamin Slough") %>%
  group_by(WS_Name, Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

# Currant Creek
# set site codes to the same site
trib_fac_sites_intersect_downstream$MonitoringLocationIdentifier[trib_fac_sites_intersect_downstream$MonitoringLocationIdentifier == "UTAHDWQ_WQX-4995312"] <- 
  "UTAHDWQ_WQX-4995310"

loading_CurrantCreek <- trib_fac_sites_intersect_downstream %>%
  ungroup() %>%
  filter(WS_Name == "Currant Creek") %>%
  group_by(WS_Name,  Month) %>%
  summarise(Flow.cfs = mean(Flow.cfs, na.rm = TRUE),
            load.TN.kgd = mean(load.TN.kgd, na.rm = TRUE), 
            load.TDN.kgd = mean(load.TDN.kgd, na.rm = TRUE),
            load.TP.kgd = mean(load.TP.kgd, na.rm = TRUE),
            load.TDP.kgd = mean(load.TDP.kgd, na.rm = TRUE),
            load.TOC.kgd = mean(load.TOC.kgd, na.rm = TRUE),
            load.DOC.kgd = mean(load.DOC.kgd, na.rm = TRUE),
            Flow.acftmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ Flow.cfs * 86400 * 31/43559.9,
                                     Month %in% c(4, 6, 9, 11) ~ Flow.cfs * 86400 * 30/43559.9,
                                     Month %in% c(2) ~ Flow.cfs * 86400 * 28/43559.9),
            load.TN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TN.kgd * 28/1000),
            load.TDN.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDN.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDN.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDN.kgd * 28/1000),
            load.TP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TP.kgd * 28/1000),
            load.TDP.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TDP.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TDP.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TDP.kgd * 28/1000),
            load.TOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.TOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.TOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.TOC.kgd * 28/1000),
            load.DOC.tonmo = case_when(Month %in% c(1, 3, 5, 7, 8, 10, 12) ~ load.DOC.kgd * 31/1000,
                                     Month %in% c(4, 6, 9, 11) ~ load.DOC.kgd * 30/1000,
                                     Month %in% c(2) ~ load.DOC.kgd * 28/1000)) %>%
  distinct()

## `summarise()` has grouped output by 'WS_Name', 'Month'. You can override using the `.groups` argument.

DMR Data

# Timp SSD
loading_TimpSSD_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "TIMPANOGOS SPECIAL SERVICE DIS") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Timp SSD")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Lindon Drain
# Calculate monthly loading from energy facility
loading_LindonDrain_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "PacifiCorp Energy") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Lindon Drain")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Powell Slough
# Calculate monthly loading from Orem WWTP
loading_PowellSloughMajor_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "OREM CITY CORP") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Powell Slough Major")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Mill Race
# Calculate monthly loading from Provo WWTP
loading_MillRace_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "PROVO CITY CORP") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Mill Race")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Hobble Creek
# Calculate loading from Springville WWTP
loading_SpringCreek_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "SPRINGVILLE- CITY OF") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Spring Creek - Springville")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Dry Creek Spanish Fork
# Calculate loading from Spanish Fork WWTP
loading_DryCreekSpanishFork_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name == "SPANISH FORK CITY CORP") %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  mutate(WS_Name = "Dry Creek - Spanish Fork")

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

# Benjamin Slough
# Calculate loading from Payson and Salem. Add loads together.

loading_BenjaminSlough_dmr <- dmr_updated_wide %>%
  ungroup() %>%
  filter(NPES_Name %in% c("SALEM CITY CORP", "PAYSON CITY CORP")) %>%
  group_by(NPES_Name, Month) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = mean(TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = mean(PO4.tonmo, na.rm = TRUE),
            load.TN.tonmo = mean(TN.tonmo, na.rm = TRUE)) %>%
  distinct() %>%
  group_by(Month, NPES_Name) %>%
  summarise(Flow.acftmo = mean(Flow.acftmo, na.rm = TRUE),
            load.TP.tonmo = sum(load.TP.tonmo, na.rm = TRUE), 
            load.TDP.tonmo = sum(load.TDP.tonmo, na.rm = TRUE),
            load.TN.tonmo = sum(load.TN.tonmo, na.rm = TRUE)) %>%
  mutate(WS_Name = "Benjamin Slough") %>%
  select(NPES_Name, Month, Flow.acftmo, load.TP.tonmo, load.TDP.tonmo, load.TN.tonmo, WS_Name)

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'Month'. You can override using the `.groups` argument.

loading_AllWatersheds_dmr <- 
  rbind(loading_BenjaminSlough_dmr, loading_DryCreekSpanishFork_dmr, 
        loading_SpringCreek_dmr, loading_MillRace_dmr,
        loading_PowellSloughMajor_dmr, loading_TimpSSD_dmr)

loading_AllWatersheds_dmr$WS_Name <- 
  factor(loading_AllWatersheds_dmr$WS_Name, 
         levels = c("Timp SSD", "Powell Slough Major",  
                    "Mill Race", "Spring Creek - Springville", 
                    "Dry Creek - Spanish Fork", "Benjamin Slough"))
loading_AllWatersheds_dmr <- loading_AllWatersheds_dmr %>%
  mutate_if(is.numeric, round, digits = 2)

## `mutate_if()` ignored the following grouping variables:
## Column `Month`

loading_AllWatersheds_dmr_annual <- loading_AllWatersheds_dmr %>%  
  group_by(NPES_Name, WS_Name) %>%
  summarise(Flow.acftyr = sum(Flow.acftmo, na.rm = TRUE),
            load.TP.tonyr = sum(load.TP.tonmo, na.rm = TRUE), 
            load.TDP.tonyr = sum(load.TDP.tonmo, na.rm = TRUE), 
            load.TN.tonyr = sum(load.TN.tonmo, na.rm = TRUE))

## `summarise()` has grouped output by 'NPES_Name'. You can override using the `.groups` argument.

loading_AllWatersheds_dmr2_annual <- loading_AllWatersheds_dmr_annual %>%
  ungroup() %>%
  group_by(WS_Name) %>%
  summarise(Flow.acftyr = sum(Flow.acftyr), 
            load.TP.tonyr = sum(load.TP.tonyr),
            load.TDP.tonyr = sum(load.TDP.tonyr),
            load.TN.tonyr = sum(load.TN.tonyr)) %>%
  mutate(OrganizationIdentifier = "DMR", 
         Flow.cfs = NA, 
         load.TDN.tonyr = NA,
         load.TOC.tonyr = NA,
         load.DOC.tonyr = NA) %>%
  select(WS_Name, Flow.cfs, load.TN.tonyr, load.TDN.tonyr, load.TP.tonyr, 
         Flow.acftyr, load.TDP.tonyr, load.TOC.tonyr, load.DOC.tonyr)

# write.csv(loading_AllWatersheds_dmr,
#           file = "./Data/Processed/loading_AllWatersheds_dmr_updated.csv",
#        row.names = FALSE)
# 
# write.csv(loading_AllWatersheds_dmr2_annual,
#           file = "./Data/Processed/loading_AllWatersheds_dmr_annual_updated.csv",
#        row.names = FALSE)

All watersheds, export

loading_AllWatersheds <- rbind(
  loading_TickvilleWash, loading_DryCreekSaratoga, 
  loading_LehiSpringCreek, loading_AmericanForkRiver, loading_TimpSSD, 
  loading_LindonDrain, loading_PowellSloughMajor, loading_ProvoRiver, 
  loading_MillRace, loading_SpringCreek, loading_HobbleCreek, loading_DryCreekSpanishFork, 
  loading_SpanishFork, loading_4000SouthDrain, loading_BenjaminSlough, 
  loading_CurrantCreek)

# calculate difference between UTAHDWQ and WFWQC load estimates
# loading_AllWatersheds_differences <- loading_AllWatersheds %>%
#   filter(WS_Name %in% c("Lehi Spring Creek", "American Fork River", "Timp SSD",
#                         "Lindon Drain", "Powell Slough Major", "Provo River", 
#                         "Mill Race", "Spring Creek - Springville", "Hobble Creek", 
#                         "Dry Creek - Spanish Fork", "Spanish Fork River", 
#                         "4000 South Drain Spanish Fork", "Benjamin Slough")) %>%
#   select(-(load.TN.kgd:load.DOC.kgd)) %>%
#   group_by(WS_Name, Month) %>%
#   summarise(load.TN.difference = load.TN.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.TN.tonmo[OrganizationIdentifier == "WFWQC_UT"], 
#             load.TDN.difference = load.TDN.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.TDN.tonmo[OrganizationIdentifier == "WFWQC_UT"], 
#             load.TP.difference = load.TP.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.TP.tonmo[OrganizationIdentifier == "WFWQC_UT"], 
#             load.TDP.difference = load.TDP.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.TN.tonmo[OrganizationIdentifier == "WFWQC_UT"], 
#             load.TOC.difference = load.TOC.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.TOC.tonmo[OrganizationIdentifier == "WFWQC_UT"], 
#             load.DOC.difference = load.DOC.tonmo[OrganizationIdentifier == "UTAHDWQ_WQX"] - 
#               load.DOC.tonmo[OrganizationIdentifier == "WFWQC_UT"])

loading_AllWatersheds_annual <- loading_AllWatersheds %>%
  group_by(WS_Name) %>%
  summarise(Flow.cfs = sum(Flow.cfs), 
            load.TN.tonyr = sum(load.TN.tonmo, na.rm = TRUE), 
            load.TDN.tonyr = sum(load.TDN.tonmo, na.rm = TRUE), 
            load.TP.tonyr = sum(load.TP.tonmo, na.rm = TRUE), 
            Flow.acftyr = sum(Flow.acftmo, na.rm = TRUE),
            load.TDP.tonyr = sum(load.TDP.tonmo, na.rm = TRUE), 
            load.TOC.tonyr = sum(load.TOC.tonmo, na.rm = TRUE), 
            load.DOC.tonyr = sum(load.DOC.tonmo, na.rm = TRUE))

loading_Total_annual <- loading_AllWatersheds_annual %>%
  ungroup() %>%
  summarise(Flow.acftyr = sum(Flow.acftyr, na.rm = TRUE), 
            load.TN.tonyr = sum(load.TN.tonyr, na.rm = TRUE), 
            load.TDN.tonyr = sum(load.TDN.tonyr, na.rm = TRUE), 
            load.TP.tonyr = sum(load.TP.tonyr, na.rm = TRUE), 
            load.TDP.tonyr = sum(load.TDP.tonyr, na.rm = TRUE), 
            load.TOC.tonyr = sum(load.TOC.tonyr, na.rm = TRUE), 
            load.DOC.tonyr = sum(load.DOC.tonyr, na.rm = TRUE))
  

# write.csv(loading_AllWatersheds, file = "./Data/Processed/loading_AllWatersheds.csv",
#        row.names = FALSE)
# write.csv(loading_AllWatersheds_annual, file = "./Data/Processed/loading_AllWatersheds_annual.csv",
#        row.names = FALSE)

# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.TN.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "TN load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")
# 
# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.TDN.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "TDN load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")
# 
# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.TP.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "TP load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")
# 
# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.TDP.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "TDP load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")
# 
# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.TOC.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "TOC load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")
# 
# ggplot(loading_AllWatersheds_differences, aes(x = Month, y = load.DOC.difference)) +
#   geom_point(alpha = 0.7, size = 2) +
#   facet_wrap(vars(WS_Name), ncol = 3, scales = "free_y") +
#   scale_x_continuous(n.breaks = 12) +
#   labs(y = "DOC load difference (UTAHDWQ-WFWQC, ton/mo)", color = "Organization") +
#   scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.7, direction = -1) +
#   theme(legend.position = "top")

# ggplot(loading_AllWatersheds, 
#        aes(x = WS_Name, y = load.TN.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   facet_wrap(vars(OrganizationIdentifier)) +
#   labs(x = "", y = "TN loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds_dmr, 
#        aes(x = WS_Name, y = load.TKN.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "TN loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(subset(loading_AllWatersheds, OrganizationIdentifier == "UTAHDWQ_WQX"),
#        aes(x = WS_Name, y = load.TDN.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "TDN loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds, 
#        aes(x = WS_Name, y = load.TP.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   facet_wrap(vars(OrganizationIdentifier)) +
#   labs(x = "", y = "TP loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds_dmr, 
#        aes(x = WS_Name, y = load.TP.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "TN loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds, 
#        aes(x = WS_Name, y = load.TDP.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   facet_wrap(vars(OrganizationIdentifier)) +
#   labs(x = "", y = "TDP loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds_dmr, 
#        aes(x = WS_Name, y = load.TDP.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "TN loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(subset(loading_AllWatersheds, OrganizationIdentifier == "UTAHDWQ_WQX"),
#        aes(x = WS_Name, y = load.TOC.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "TOC loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(subset(loading_AllWatersheds, OrganizationIdentifier == "UTAHDWQ_WQX"),
#        aes(x = WS_Name, y = load.DOC.tonmo, fill = as.factor(Month))) +
#   geom_bar(stat = "identity") +
#     scale_fill_viridis_d(option = "magma") +
#   labs(x = "", y = "DOC loading (tons/year)", fill = "Month") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds_annual_combined, 
#        aes(x = WS_Name, y = load.TP.tonyr, fill = OrganizationIdentifier)) +
#   geom_bar(stat = "identity", position = position_dodge(preserve = "single")) +
#   scale_fill_viridis_d(option = "magma", begin = 0.2, end = 0.8) +
#   labs(x = "", y = "TP loading (tons/year)", fill = "Source") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
# 
# ggplot(loading_AllWatersheds_annual_combined, 
#        aes(x = WS_Name, y = load.TN.tonyr, fill = OrganizationIdentifier)) +
#   geom_bar(stat = "identity", position = position_dodge(preserve = "single")) +
#   scale_fill_viridis_d(option = "magma", begin = 0.2, end = 0.8) +
#   labs(x = "", y = "TN loading (tons/year)", fill = "Source") +
#   theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
#

C-Q

Precipitation and evaporation

Appendix

Nutrient concentrations and flow across watersheds

Watersheds w/ UDWQ only, no WWTP

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Tickville Wash"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Tickville Wash") +
  scale_fill_manual(values = c("#f27162")) +
  theme(legend.position = "top")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Tickville.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Saratoga"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Dry Creek - Saratoga") +
  scale_fill_manual(values = c("#f27162")) +
  theme(legend.position = "top")

## Warning: Removed 26 rows containing non-finite values (stat_boxplot).

## Warning: Removed 26 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_DrySaratoga.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 26 rows containing non-finite values (stat_boxplot).

## Warning: Removed 26 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Currant Creek"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Currant Creek") +
  scale_fill_manual(values = c("#f27162", "#db4a69")) +
  theme(legend.position = "top")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Currant.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Tickville Wash"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Tickville Wash") +
  scale_fill_manual(values = c("#f27162")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Tickville.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 15 rows containing non-finite values (stat_boxplot).

## Warning: Removed 15 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Saratoga"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Dry Creek - Saratoga") +
  scale_fill_manual(values = c("#f27162")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 26 rows containing non-finite values (stat_boxplot).

## Warning: Removed 26 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_DrySaratoga.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 26 rows containing non-finite values (stat_boxplot).

## Warning: Removed 26 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Currant Creek"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Currant Creek") +
  scale_fill_manual(values = c("#f27162", "#db4a69")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 16 rows containing non-finite values (stat_boxplot).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Currant.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 16 rows containing non-finite values (stat_boxplot).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

Flow

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Tickville Wash"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Tickville Wash") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Tickville.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Saratoga"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Dry Creek - Saratoga") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_DrySaratoga.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Currant Creek"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Currant Creek") +
  scale_fill_manual(values = c("#f27162", "#db4a69")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 14 rows containing non-finite values (stat_boxplot).

## Warning: Removed 14 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Currant.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 14 rows containing non-finite values (stat_boxplot).

## Warning: Removed 14 rows containing non-finite values (stat_summary).

Watersheds w/ UDWQ and WFWQC, no WWTP

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lehi Spring Creek"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Lehi Spring Creek") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 52 rows containing non-finite values (stat_boxplot).

## Warning: Removed 52 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Lehi.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 52 rows containing non-finite values (stat_boxplot).

## Warning: Removed 52 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "American Fork River"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "American Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 59 rows containing non-finite values (stat_boxplot).

## Warning: Removed 59 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_American.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 59 rows containing non-finite values (stat_boxplot).

## Warning: Removed 59 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lindon Drain"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Lindon Drain") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 55 rows containing non-finite values (stat_boxplot).

## Warning: Removed 55 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Lindon.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 55 rows containing non-finite values (stat_boxplot).

## Warning: Removed 55 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Provo River"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Provo River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 57 rows containing non-finite values (stat_boxplot).

## Warning: Removed 57 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Provo.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 57 rows containing non-finite values (stat_boxplot).

## Warning: Removed 57 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Hobble Creek"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Hobble Creek") +
  scale_fill_manual(values = c("#fc9d6f",  "#8a2880", "#3b1165")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 63 rows containing non-finite values (stat_boxplot).

## Warning: Removed 63 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Hobble.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 63 rows containing non-finite values (stat_boxplot).

## Warning: Removed 63 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spanish Fork River"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  ylim(0, 2000) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Spanish Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 57 rows containing non-finite values (stat_boxplot).

## Warning: Removed 57 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Spanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 57 rows containing non-finite values (stat_boxplot).

## Warning: Removed 57 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "4000 South Drain Spanish Fork"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "4000 South Drain Spanish Fork") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 56 rows containing non-finite values (stat_boxplot).

## Warning: Removed 56 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_4000.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 56 rows containing non-finite values (stat_boxplot).

## Warning: Removed 56 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lehi Spring Creek"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Lehi Spring Creek") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Lehi.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "American Fork River"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "American Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 56 rows containing non-finite values (stat_boxplot).

## Warning: Removed 56 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_American.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 56 rows containing non-finite values (stat_boxplot).

## Warning: Removed 56 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lindon Drain"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Lindon Drain") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Lindon.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Provo River"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Provo River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Provo.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Hobble Creek"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Hobble Creek") +
  scale_fill_manual(values = c("#fc9d6f",  "#8a2880", "#3b1165")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")+
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 55 rows containing non-finite values (stat_boxplot).

## Warning: Removed 55 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Hobble.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 55 rows containing non-finite values (stat_boxplot).

## Warning: Removed 55 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spanish Fork River"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Spanish Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 33 rows containing non-finite values (stat_boxplot).

## Warning: Removed 33 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Spanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 33 rows containing non-finite values (stat_boxplot).

## Warning: Removed 33 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "4000 South Drain Spanish Fork"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "4000 South Drain Spanish Fork") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_4000.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

Flow

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lehi Spring Creek"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Lehi Spring Creek") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Lehi.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "American Fork River"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "American Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 34 rows containing non-finite values (stat_boxplot).

## Warning: Removed 34 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_American.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 34 rows containing non-finite values (stat_boxplot).

## Warning: Removed 34 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lindon Drain"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Lindon Drain") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Lindon.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Provo River"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Provo River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Provo.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Hobble Creek"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Hobble Creek") +
  scale_fill_manual(values = c("#fc9d6f",  "#8a2880", "#3b1165")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

## Warning: Removed 1 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Hobble.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

## Warning: Removed 1 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spanish Fork River"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Spanish Fork River") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Spanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "4000 South Drain Spanish Fork"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "4000 South Drain Spanish Fork") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 11 rows containing non-finite values (stat_boxplot).

## Warning: Removed 11 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_4000.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 11 rows containing non-finite values (stat_boxplot).

## Warning: Removed 11 rows containing non-finite values (stat_summary).

Watersheds w/ UDWQ and WFWQC, plus WWTP

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Timp SSD"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Timp SSD") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 68 rows containing non-finite values (stat_boxplot).

## Warning: Removed 68 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_TimpSSD.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 68 rows containing non-finite values (stat_boxplot).

## Warning: Removed 68 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Powell Slough Major"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Powell Slough Major") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Powell.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 49 rows containing non-finite values (stat_boxplot).

## Warning: Removed 49 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Mill Race"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Mill Race") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#8a2880", "#3b1165")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 66 rows containing non-finite values (stat_boxplot).

## Warning: Removed 66 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_MillRace.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 66 rows containing non-finite values (stat_boxplot).

## Warning: Removed 66 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spring Creek - Springville"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Spring Creek - Springville") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 32 rows containing non-finite values (stat_boxplot).

## Warning: Removed 32 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Spring.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 32 rows containing non-finite values (stat_boxplot).

## Warning: Removed 32 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Spanish Fork"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Dry Creek - Spanish Fork") +
  scale_fill_manual(values = c("#fc9d6f", "#f27162", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 3 rows containing non-finite values (stat_boxplot).

## Warning: Removed 3 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_DrySpanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 3 rows containing non-finite values (stat_boxplot).

## Warning: Removed 3 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Benjamin Slough"),
       aes(x = as.factor(Month), y = TN, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 200, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 700, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TN ("*mu*"g/L)"), fill = "Site", title = "Benjamin Slough") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 48 rows containing non-finite values (stat_boxplot).

## Warning: Removed 48 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TNconc_Benjamin.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 48 rows containing non-finite values (stat_boxplot).

## Warning: Removed 48 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Timp SSD"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Timp SSD") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 38 rows containing non-finite values (stat_boxplot).

## Warning: Removed 38 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_TimpSSD.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 38 rows containing non-finite values (stat_boxplot).

## Warning: Removed 38 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Powell Slough Major"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Powell Slough Major") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 21 rows containing non-finite values (stat_boxplot).

## Warning: Removed 21 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Powell.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 21 rows containing non-finite values (stat_boxplot).

## Warning: Removed 21 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Mill Race"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Mill Race") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#8a2880", "#3b1165")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 29 rows containing non-finite values (stat_boxplot).

## Warning: Removed 29 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_MillRace.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 29 rows containing non-finite values (stat_boxplot).

## Warning: Removed 29 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spring Creek - Springville"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Spring Creek - Springville") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 8 rows containing non-finite values (stat_boxplot).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Spring.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 8 rows containing non-finite values (stat_boxplot).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Spanish Fork"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Dry Creek - Spanish Fork") +
  scale_fill_manual(values = c("#fc9d6f", "#f27162", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

## Warning: Removed 2 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_DrySpanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

## Warning: Removed 2 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Benjamin Slough"),
       aes(x = as.factor(Month), y = TP, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  geom_hline(yintercept = 3, lty = 2, color = "#f27162") +
  geom_hline(yintercept = 21, lty = 2, color = "#62197d") +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = expression("TP ("*mu*"g/L)"), fill = "Site", title = "Benjamin Slough") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

## Warning: Removed 23 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/TPconc_Benjamin.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

## Warning: Removed 23 rows containing non-finite values (stat_summary).

Flow

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Lehi Spring Creek"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier, 
           color = Backwater)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Lehi Spring Creek") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  scale_color_manual(values = c("black", "gray"), guide = "none") +
  theme(legend.position = "top")

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Lehi.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 18 rows containing non-finite values (stat_boxplot).

## Warning: Removed 18 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Timp SSD"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10(limits = c(10, 100)) +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Timp SSD") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 52 rows containing non-finite values (stat_boxplot).

## Warning: Removed 52 rows containing non-finite values (stat_summary).

## Warning: Removed 1 rows containing missing values (geom_text).

ggsave("./Code/Output/Flow_TimpSSD.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 52 rows containing non-finite values (stat_boxplot).

## Warning: Removed 52 rows containing non-finite values (stat_summary).

## Warning: Removed 1 rows containing missing values (geom_text).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Powell Slough Major"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Powell Slough Major") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 8 rows containing non-finite values (stat_boxplot).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Powell.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 8 rows containing non-finite values (stat_boxplot).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Mill Race"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Mill Race") +
  scale_fill_manual(values = c("#fc9d6f", "#db4a69", "#8a2880", "#3b1165")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 13 rows containing non-finite values (stat_boxplot).

## Warning: Removed 13 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_MillRace.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 13 rows containing non-finite values (stat_boxplot).

## Warning: Removed 13 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Spring Creek - Springville"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Spring Creek - Springville") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

## Warning: Removed 2 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Spring.jpg", height = 4, width = 5, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

## Warning: Removed 2 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Dry Creek - Spanish Fork"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Dry Creek - Spanish Fork") +
  scale_fill_manual(values = c("#fc9d6f", "#f27162", "#db4a69", "#62197d")) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(nrow = 2))

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_DrySpanish.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing non-finite values (stat_summary).

ggplot(subset(trib_fac_sites_intersect_downstream, WS_Name == "Benjamin Slough"),
       aes(x = as.factor(Month), y = Flow.cfs, fill = MonitoringLocationIdentifier)) +
  geom_boxplot(position = position_dodge(preserve = "single"), alpha = 0.8) +
  scale_y_log10() +
  stat_summary(fun.data = give.n, geom = "text", size = 3,
                  position = position_dodge(width = 0.75)) +
  labs(x = "Month", y = "Flow (cfs)", fill = "Site", title = "Benjamin Slough") +
  scale_fill_manual(values = c("#f27162", "#62197d")) +
  theme(legend.position = "top")

## Warning: Removed 13 rows containing non-finite values (stat_boxplot).

## Warning: Removed 13 rows containing non-finite values (stat_summary).

ggsave("./Code/Output/Flow_Benjamin.jpg", height = 4, width = 5, units = "in")

## Warning: Removed 13 rows containing non-finite values (stat_boxplot).

## Warning: Removed 13 rows containing non-finite values (stat_summary).

UtahLake_TributaryData

Kateri Salk

Details

Load data

Data processing

Watersheds

Tributary data

Find non-sampled dates

Facility data

Combine trib and fac data

USGS gage data for missing discharge

Create spatial datasets for tributary and facility datasets

Jordan River Outflow

Groundwater (spring and well) data

Precipitation and Evapotranspiration data

DMR data

Visualization

Maps

Water Quality Time Series

Flow time series

Load time series

DMR Flow, concentration, and load

Monthly loads by watershed

Tickville Wash

Dry Creek - Saratoga

Lehi Spring Creek

American Fork River

Timp SSD

Lindon Drain

Powell Slough Major

Provo River

Mill Race

Spring Creek - Springville

Hobble Creek

Dry Creek - Spanish Fork

Spanish Fork River

4000 South Drain Spanish Fork

Benjamin Slough

Currant Creek

Monthly and annual load across watersheds

Monitored tributaries

DMR Data

All watersheds, export

C-Q

Precipitation and evaporation

Appendix

Nutrient concentrations and flow across watersheds

Watersheds w/ UDWQ only, no WWTP

Watersheds w/ UDWQ and WFWQC, no WWTP

Watersheds w/ UDWQ and WFWQC, plus WWTP

Tickville Wash

Dry Creek - Saratoga

Lehi Spring Creek

American Fork River

Timp SSD

Lindon Drain

Powell Slough Major

Provo River

Mill Race

Spring Creek - Springville

Hobble Creek

Dry Creek - Spanish Fork

Spanish Fork River

4000 South Drain Spanish Fork

Benjamin Slough

Currant Creek

Spare code